### Install and Initialize QuantumToolbox.jl

Source: https://context7.com/qutip/quantumtoolbox.jl/llms.txt

Demonstrates how to install the package using Julia's Pkg manager and verify the installation version.

```julia
using Pkg
Pkg.add("QuantumToolbox")

# Load the package
using QuantumToolbox

# Check version information
QuantumToolbox.versioninfo()
```

--------------------------------

### Install and Initialize QuantumToolbox.jl

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/README.md

Instructions for installing the package via the Julia Pkg manager and verifying the installation using built-in version information tools.

```julia
using Pkg
Pkg.add("QuantumToolbox")
```

```julia-repl
(1.10) pkg> add QuantumToolbox
```

```julia
using QuantumToolbox
QuantumToolbox.versioninfo()
QuantumToolbox.about()
```

--------------------------------

### Start Local Vitepress Site with Make

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/README.md

Starts a local Vitepress development server for the QuantumToolbox.jl documentation using the `make vitepress` command. This allows for live previewing of documentation changes. Requires Node.js and npm.

```shell
make vitepress
```

--------------------------------

### Start Local Vitepress Site Manually

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/README.md

Starts a local Vitepress development server for the QuantumToolbox.jl documentation by running the `docs:dev` script using npm. This allows for live previewing of documentation changes. Requires Node.js and npm.

```shell
npm --prefix docs run docs:dev
```

--------------------------------

### Install QuantumToolbox.jl Package

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/index.md

Installs the QuantumToolbox.jl package using Julia's Pkg manager. This can be done directly in the Julia REPL or Pkg REPL. Requires Julia 1.10+.

```julia
using Pkg
Pkg.add("QuantumToolbox")
```

```julia-repl
(1.10) pkg> add QuantumToolbox
```

--------------------------------

### Setup Makie for QuantumToolbox.jl Visualization

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Initializes Makie.jl for plotting within the QuantumToolbox.jl environment. It enables SVG output and imports necessary libraries for visualization.

```julia
using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
using QuantumToolbox
```

--------------------------------

### QuantumToolbox.jl: Harmonic Oscillator Steady State Example

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/steadystate.md

A comprehensive example demonstrating the `steadystate` function for a harmonic oscillator in a thermal bath. It includes setting up the system, defining collapse operators, calculating the steady-state solution, and comparing it with solutions from Monte-Carlo and Master Equation solvers. Plots are generated using CairoMakie.

```julia
using QuantumToolbox
using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())

# Define parameters
N = 20  # number of basis states to consider
a = destroy(N)
H = a' * a
ψ0 = basis(N, 10)  # initial state
κ = 0.1  # coupling to oscillator
n_th = 2 # temperature with average of 2 excitations

# collapse operators 
c_ops = [
    sqrt(κ * (n_th + 1)) * a, # emission
    sqrt(κ *  n_th     ) * a' # absorption
]

# find steady-state solution
ρ_ss = steadystate(H, c_ops)

# find expectation value for particle number in steady state
e_ops = [a' * a]
exp_ss = real(expect(e_ops[1], ρ_ss))

tlist = LinRange(0, 50, 100)

# monte-carlo
sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=100, progress_bar=false)
exp_mc = real(sol_mc.expect[1, :])

# master eq.
sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=false)
exp_me = real(sol_me.expect[1, :])

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], 
    title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$",
    xlabel = "Time", 
    ylabel = "Number of excitations",
)
lines!(ax, tlist, exp_mc, label = "Monte-Carlo", linewidth = 2, color = :blue)
lines!(ax, tlist, exp_me, label = "Master Equation", linewidth = 2, color = :orange, linestyle = :dash)
lines!(ax, tlist, exp_ss .* ones(length(tlist)), label = "Steady State", linewidth = 2, color = :red)
axislegend(ax, position = :rt)

fig
```

--------------------------------

### QuantumToolbox Setup for Time Evolution

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/solution.md

Sets up the QuantumToolbox environment for time evolution simulations, including necessary using declarations and enabling SVG output for plotting.

```julia
using QuantumToolbox

using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())

```

--------------------------------

### Build Documentation with Make

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/README.md

Builds the documentation for the QuantumToolbox.jl project using the `make docs` command. This method requires Node.js and npm to be installed. It generates documentation files in the `docs/build/1/` directory.

```shell
make docs
```

--------------------------------

### Full Example: Time-Dependent Hamiltonian and Collapse Operators (Julia)

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/time_dependent.md

Sets up and solves a quantum system with a time-dependent Hamiltonian and time-dependent collapse operators using `mesolve` and `mcsolve`. This example demonstrates defining operators, coefficient functions, initial states, time lists, and observables for a multi-level system.

```julia
N = 2 # Set where to truncate Fock state for cavity

# basis states for Atom
u = basis(3, 0) #       u state
e = basis(3, 1) # excited state
g = basis(3, 2) #  ground state

# operators
g_e = tensor(qeye(N), g * e')  # |g><e|
u_e = tensor(qeye(N), u * e')  # |u><e|
u_u = tensor(qeye(N), u * u')  # |u><u|
g_g = tensor(qeye(N), g * g')  # |g><g|
a = tensor(destroy(N), qeye(3))

# Hamiltonian
g = 5  # coupling strength
H0 = -g * (g_e' * a + a' * g_e)
H1 = u_e' + u_e
f1(p, t) = 9 * exp(-(t / 5)^2)
H_t = QobjEvo(
    (
        H0,
        (H1, f1)
    )
)

# Build collapse operators
κ = 1.5 # Cavity decay rate
γ0 = 6   # Atomic decay rate
c_ops = [
    sqrt(κ) * a,
    sqrt(4 * γ0 / 9) * g_e, # Use Rb branching ratio of 4/9 e -> g
    sqrt(5 * γ0 / 9) * u_e  # 5/9 e -> u
]

tlist = LinRange(0, 4, 200) # Define time vector
ψ0 = tensor(basis(N, 0), u) # Define initial state

# Build observables
e_ops = [
    a' * a,
    u_u,
    g_g
]

# solve dynamics 
exp_me = mesolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, progress_bar = Val(false)).expect
exp_mc = mcsolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, ntraj = 100, progress_bar = Val(false)).expect
```

--------------------------------

### Generate Quantum States

Source: https://context7.com/qutip/quantumtoolbox.jl/llms.txt

Provides examples for creating various quantum states including Fock states, coherent states, density matrices, and spin states.

```julia
using QuantumToolbox

# Fock state (number state) |n>
ψ0 = fock(5, 0)
ψ = basis(5, 2)

# Coherent state |α>
α = 0.5 - 0.5im
ψ_coh = coherent(10, α)

# Density matrices
ρ_fock = fock_dm(5, 2)
ρ_thermal = thermal_dm(10, 2.0)

# Spin states
spin_up = basis(2, 0)
bell = bell_state(0, 0)
```

--------------------------------

### Build Documentation Manually

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/README.md

Manually builds the documentation for the QuantumToolbox.jl project by executing the `docs/make.jl` script using Julia. This method requires Node.js and npm to be installed. It generates documentation files in the `docs/build/1/` directory.

```julia
julia --project=docs docs/make.jl
```

--------------------------------

### Setup for Master Equation Solver (Julia)

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/mesolve.md

This snippet sets up the necessary environment for using the master equation solver in QuantumToolbox.jl. It imports the required library and configures plotting for SVG output.

```julia
using QuantumToolbox

using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
```

--------------------------------

### Initialize Quantum Objects in Julia and Python

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/QuantumObject/QuantumObject.md

Demonstrates how to instantiate QuantumObject (Julia) or Qobj (Python) using various array inputs. It includes examples of defining Hilbert space dimensions using Tuples or StaticArrays for improved performance.

```julia
using QuantumToolbox
QuantumObject([1, 2, 3, 4, 5])
M = rand(ComplexF64, 4, 4)
Qobj(M, dims = (2, 2))
Qobj(rand(4, 4), type = SuperOperator())
```

```python
from qutip import Qobj
Qobj([1, 2, 3, 4, 5])
Qobj([1, 2, 3, 4, 5])
Qobj(rand(4, 4))
```

--------------------------------

### Master Equation with Density Matrix Initial State (Julia)

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/mesolve.md

Illustrates using `mesolve` when the initial state is provided as a density matrix. This example also shows how to specify the time points for saving the solution and calculating expectation values.

```julia
tlist = [0, 5, 10]
state0 = ket2dm(basis(2, 0))  # density matrix
sol = mesolve(H, state0, tlist, e_ops = [sigmay()], saveat = tlist)
```

--------------------------------

### Running Monte Carlo simulation with full trajectory storage

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/solution.md

Example of executing mcsolve with keep_runs_results=Val(true). This configuration stores the full history for every trajectory.

```julia
sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(true), progress_bar=Val(false))
size(sol_mc2.expect)
```

--------------------------------

### QuantumToolbox Time Evolution Solution Example (Julia)

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/solution.md

Demonstrates simulating quantum system dynamics using QuantumToolbox.jl's sesolve function. It defines a Hamiltonian, initial state, expectation operators, and a time list to generate a time evolution solution.

```julia
H = 0.5 * sigmay()
ψ0 = basis(2, 0)
e_ops = [
    proj(basis(2, 0)),
    proj(basis(2, 1)),
    basis(2, 0) * basis(2, 1)'
]
tlist = LinRange(0, 10, 100)
sol = sesolve(H, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false))

```

--------------------------------

### Constructing Liouvillians with QobjEvo and Lindblad Dissipators

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/time_dependent.md

This example shows how to construct a time-dependent Liouvillian superoperator using QobjEvo. It demonstrates both direct construction using spre and spost, and the optimized approach using the liouvillian function.

```julia
using QuantumOptics

coef(p, t) = sin(π * t)
Ht = QobjEvo(sigmaz(), coef)

coef2(p, t) = exp(-t)
c_op = QobjEvo(destroy(2), coef2)

L1 = -1im * (spre(Ht) - spost(Ht'))
L1 += lindblad_dissipator(c_op)

# Or equivalently:
L2 = liouvillian(Ht, [c_op])

t = rand()
println(L1(t) == L2(t))
```

--------------------------------

### Load and Check QuantumToolbox.jl Version

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/index.md

Loads the QuantumToolbox.jl package into the current Julia session and displays version information or general package details. This is useful for verifying installation and understanding the package's current state.

```julia
using QuantumToolbox
QuantumToolbox.versioninfo()
QuantumToolbox.about()
```

--------------------------------

### Simulate Homodyne Detection with ssesolve in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/stochastic.md

This example demonstrates how to set up and solve the stochastic Schrödinger equation for an optical cavity. It defines the Hamiltonian, initial coherent state, and stochastic collapse operators before executing the simulation with ssesolve.

```julia
using QuantumToolbox

# parameters
N = 20         # Fock space dimension
Δ = 5 * 2 * π  # cavity detuning
κ = 2          # cavity decay rate
α = 4          # intensity of initial state
ntraj = 500    # number of trajectories

tlist = 0:0.0025:1

# operators
a = destroy(N)
x = a + a'
H = Δ * a' * a

# initial state
ψ0 = coherent(N, √α)

# temperature with average of 0 excitations (absolute zero)
n_th = 0
sc_ops = [√(κ * (n_th + 1)) * a]

# Solve using ssesolve
sse_sol = ssesolve(
    H,
    ψ0,
    tlist,
    sc_ops,
    e_ops = [x],
    ntraj = ntraj,
    store_measurement = Val(true),
)

measurement_avg = sum(sse_sol.measurement, dims=2) / size(sse_sol.measurement, 2)
measurement_avg = dropdims(measurement_avg, dims=2)
```

--------------------------------

### CPU Computation: Quantum Harmonic Oscillator Simulation with QuantumToolbox.jl

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/getting_started/brief_example.md

This snippet demonstrates CPU computation for a quantum harmonic oscillator using QuantumToolbox.jl. It defines the Hamiltonian, sets up the Lindblad master equation with damping, and computes the time evolution of the system's density matrix and expectation values.

```julia
using QuantumToolbox

N = 20 # cutoff of the Hilbert space dimension
ω = 1.0 # frequency of the harmonic oscillator

a = destroy(N) # annihilation operator

H = ω * a' * a

γ = 0.1 # damping rate

ψ0 = fock(N, 3) # initial state

tlist = range(0, 10, 100) # time list

c_ops = [sqrt(γ) * a]
e_ops = [a' * a]

sol = mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops)
```

--------------------------------

### GPU Computation: Quantum Harmonic Oscillator Simulation with QuantumToolbox.jl and CUDA.jl

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/getting_started/brief_example.md

This snippet shows how to perform the same quantum harmonic oscillator simulation on the GPU using QuantumToolbox.jl and CUDA.jl. The primary difference is the use of the `cu()` function to transfer quantum objects to the GPU.

```julia
using QuantumToolbox
using CUDA
CUDA.allowscalar(false) # Avoid unexpected scalar indexing

# Assuming N, ω, γ, tlist are defined as in the CPU example

a_gpu = cu(destroy(N)) # The only difference in the code is the cu() function

H_gpu = ω * a_gpu' * a_gpu

ψ0_gpu = cu(fock(N, 3))

c_ops = [sqrt(γ) * a_gpu]
e_ops = [a_gpu' * a_gpu]

sol = mesolve(H_gpu, ψ0_gpu, tlist, c_ops, e_ops = e_ops)
```

--------------------------------

### Lindblad Master Equation Solver (mesolve) in Julia

Source: https://context7.com/qutip/quantumtoolbox.jl/llms.txt

Explains how to use `mesolve` in QuantumToolbox.jl to simulate open quantum systems evolving under the Lindblad master equation, including dissipation and decoherence. It shows examples with both ket and density matrix initial states.

```julia
using QuantumToolbox

# Harmonic oscillator with damping
N = 20
a = destroy(N)
H = a' * a                    # Hamiltonian
ψ0 = fock(N, 10)              # initial Fock state |10>

# Thermal bath parameters
κ = 0.1                       # coupling rate
n_th = 2                      # thermal occupation

# Collapse operators for thermal bath
c_ops = [
    sqrt(κ * (n_th + 1)) * a,  # emission (cooling)
    sqrt(κ * n_th) * a'        # absorption (heating)
]

tlist = range(0, 50, 100)
e_ops = [a' * a]              # number operator

# Solve master equation
sol = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
photon_number = real(sol.expect[1, :])

# With density matrix initial state
ρ0 = ket2dm(ψ0)
sol = mesolve(H, ρ0, tlist, c_ops, e_ops=e_ops)

# Generate Liouvillian superoperator manually
L = liouvillian(H, c_ops)
sol = mesolve(L, ψ0, tlist, e_ops=e_ops)
```

--------------------------------

### QuantumToolbox.jl: Leveraging Intel MKL for LU Decomposition

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/steadystate.md

Demonstrates how to utilize Intel MKL's optimized LU solvers (Pardiso and standard LU) through `LinearSolve.jl` for potentially faster direct steady-state solutions. This requires loading `Pardiso.jl` or `MKL_jll.jl`.

```julia
using QuantumToolbox
using LinearSolve # must be loaded

using Pardiso
solver = SteadyStateLinearSolver(alg = MKLPardisoFactorize())

using MKL_jll
solver = SteadyStateLinearSolver(alg = MKLLUFactorization())
```

--------------------------------

### QuantumToolbox.jl: Specifying a Solver for steadystate

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/steadystate.md

Shows how to explicitly select a solver for the `steadystate` function using the `solver` keyword argument. This allows for customization of the numerical method used to find the steady-state solution.

```julia
ρ_ss = steadystate(H, c_ops; solver = SteadyStateLinearSolver())
```

--------------------------------

### Demonstrating Type Instability in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/getting_started/type_stability.md

An example of a type-unstable function where the return type depends on the input value, causing runtime type inference issues.

```julia
function foo(x)
    if x > 0
        return 1
    else
        return -1.0
    end
end
```

--------------------------------

### Convergence Study of Monte Carlo Solver (QuantumToolbox.jl)

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/mcsolve.md

Demonstrates the convergence of the Monte Carlo solver by simulating the same system with varying numbers of trajectories (`ntraj`). It plots the expectation values for 1, 10, and 100 trajectories to visualize how the results approach a stable solution. This helps in determining an adequate number of trajectories for accurate simulations.

```julia
using QuantumToolbox
using CairoMakie

tlist = LinRange(0.0, 10.0, 200)

ψ0 = tensor(fock(2, 0), fock(10, 8))
a  = tensor(qeye(2), destroy(10))
σm = tensor(destroy(2), qeye(10))
H = 2 * π * a' * a + 2 * π * σm' * σm + 2 * π * 0.25 * (σm * a' + σm' * a)

c_ops = [sqrt(0.1) * a]
e_ops = [a' * a]

sol_1   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 1)
sol_10  = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 10)
sol_100 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 100)

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1],
    xlabel = "Time",
    ylabel = "Expectation values",
    title = "Monte Carlo time evolution",
)
lines!(ax, tlist, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
lines!(ax, tlist, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
lines!(ax, tlist, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)

axislegend(ax, position = :rt)

fig
```

--------------------------------

### Running Monte Carlo simulation with averaged results

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/solution.md

Example of executing mcsolve with keep_runs_results=Val(false). This configuration averages results across all trajectories before storage.

```julia
tlist = LinRange(0, 1, 11)
c_ops = (destroy(2),)
e_ops = (num(2),)
sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(false), progress_bar=Val(false))
size(sol_mc1.expect)
```

--------------------------------

### QuantumToolbox.jl: Using KrylovJL Solvers with Preconditioners

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/steadystate.md

Illustrates how to configure the `SteadyStateLinearSolver` with iterative algorithms like `KrylovJL_GMRES` and apply preconditioners for potentially faster convergence, especially for large systems. Note that preconditioners are not guaranteed to be effective for non-Hermitian Liouvillians.

```julia
SteadyStateLinearSolver(alg = KrylovJL_GMRES(; precs = (A, p) -> (ilu(A, τ = 0.01), I)))
```

--------------------------------

### Generate Density Matrices for Quantum States

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Generates density matrices for coherent, thermal, and Fock states using QuantumToolbox.jl. These are used as input for subsequent visualization functions.

```julia
N = 20

ρ_coherent = coherent_dm(N, sqrt(2))
ρ_thermal = thermal_dm(N, 2)
ρ_fock = fock_dm(N, 2)
nothing # hide
```

--------------------------------

### Qobj Inner Product: Bra * Ket

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/QuantumObject/QuantumObject.md

Shows the inner product calculation (<v|u>) when a Bra vector is multiplied by a Ket vector. This is the reverse order of the outer product example.

```julia
v' * u
```

--------------------------------

### Initialize environment for QuantumToolbox.jl

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/getting_started/logo.md

Loads the necessary QuantumToolbox and CairoMakie packages and configures the plotting backend to output SVG images.

```julia
using QuantumToolbox
using CairoMakie
CairoMakie.activate!(type = "svg", pt_per_unit = 1)
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
```

--------------------------------

### Plot Fock-Basis Probability Distribution

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Visualizes the Fock-basis probability distribution of quantum states using QuantumToolbox.jl's plot_fock_distribution function and Makie.jl. It displays distributions for coherent, thermal, and Fock states.

```julia
fig = Figure(size = (600, 350))

# only show several x-labels
fock_numbers = [
    "0", "", "", "", "", "5", "", "", "", "",
    "10", "", "", "", "", "15", "", "", "", "19",
]

_, ax1, _ = plot_fock_distribution(ρ_coherent; location=fig[1,1], fock_numbers=fock_numbers)
ax1.title = "Coherent state"

_, ax2, _ =plot_fock_distribution(ρ_thermal; location=fig[1,2], fock_numbers=fock_numbers)
ax2.title = "Thermal state"

_, ax3, _ = plot_fock_distribution(ρ_fock; location=fig[1,3], fock_numbers=fock_numbers)
ax3.title = "Fock state"

fig
```

--------------------------------

### Tensor Products and Partial Traces in Julia

Source: https://context7.com/qutip/quantumtoolbox.jl/llms.txt

Illustrates how to construct composite quantum systems using tensor products and how to perform partial traces to focus on specific subsystems. This is crucial for multi-particle systems and analyzing entanglement.

```julia
using QuantumToolbox

# Tensor product of states
ψ1 = basis(2, 0)  # |0>
ψ2 = basis(2, 1)  # |1>
ψ_composite = tensor(ψ1, ψ2)  # |0>⊗|1> = |01>

# Tensor product of operators
σz_1 = tensor(sigmaz(), qeye(2))  # σz ⊗ I
σz_2 = tensor(qeye(2), sigmaz())  # I ⊗ σz

# Multiple subsystems
three_qubits = tensor(basis(2,0), basis(2,1), basis(2,0))

# Construct composite Hamiltonian (two coupled qubits)
H = tensor(sigmaz(), qeye(2)) +
    tensor(qeye(2), sigmaz()) +
    0.05 * tensor(sigmax(), sigmax())

# Jaynes-Cummings model (qubit + cavity)
N_cavity = 10
ωa, ωc, g = 1.0, 1.25, 0.75
a = tensor(qeye(2), destroy(N_cavity))   # cavity operator
σm = tensor(destroy(2), qeye(N_cavity))  # qubit lowering
σz = tensor(sigmaz(), qeye(N_cavity))
H_JC = 0.5*ωa*σz + ωc*a'*a + g*(a'*σm + a*σm')

# Partial trace - trace out subsystems
ψ = tensor(basis(2, 0), basis(2, 1), normalize(basis(2, 0) + basis(2, 1)))
ρ1 = ptrace(ψ, 1)        # keep first subsystem, trace out 2nd and 3rd
ρ13 = ptrace(ψ, (1, 3))  # keep 1st and 3rd, trace out 2nd
```

--------------------------------

### Create and Manipulate QuantumObjects (Qobj)

Source: https://context7.com/qutip/quantumtoolbox.jl/llms.txt

Shows how to instantiate Qobj structures from arrays or matrices and access their internal properties like dimensions and data.

```julia
using QuantumToolbox

# Create a ket vector from an array
ψ = Qobj([1, 2, 3, 4, 5])

# Create an operator from a matrix
M = Qobj([1 2; 3 4])

# Create operator with specific dimensions for composite systems
M = Qobj(rand(ComplexF64, 4, 4), dims=(2, 2))

# Access properties
println(M.data)
println(M.type)
println(M.dims)
println(size(M))
println(ishermitian(M))
```

--------------------------------

### Simulating Harmonic Oscillator Decay with Time-Dependent Decay Rate

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/time_dependent.md

This example simulates the decay of a harmonic oscillator with both time-independent and time-dependent decay rates using mesolve. It then plots the expectation value of the number operator for both cases.

```julia
using QuantumOptics, CairoMakie

N = 10  # number of basis states
a = destroy(N)
H = a' * a

ψ0 = basis(N, 9)  # initial state

γ0 = 0.5
γ1(p, t) = sqrt(γ0 * exp(-t))
c_ops_ti = [sqrt(γ0) * a]   # time-independent collapse term
c_ops_td = [QobjEvo(a, γ1)] # time-dependent   collapse term

e_ops = [a' * a]

tlist = LinRange(0, 10, 100)
exp_ti = mesolve(H, ψ0, tlist, c_ops_ti, e_ops = e_ops; progress_bar = Val(false)).expect[1,:]
exp_td = mesolve(H, ψ0, tlist, c_ops_td, e_ops = e_ops; progress_bar = Val(false)).expect[1,:]

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = L"Time $t$", ylabel = L"$\langle \hat{a}^\dagger \hat{a} \rangle$")
lines!(ax, tlist, real(exp_ti), label = L"$\gamma_0$", linestyle = :solid)
lines!(ax, tlist, real(exp_td), label = L"$\gamma_1(t) = \gamma_0 e^{-t}$", linestyle = :dash)

axislegend(ax, position = :rt)

fig
```

--------------------------------

### Git Commands for Contribution to QuantumToolbox.jl

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/resources/contributing.md

This snippet demonstrates essential Git commands for contributing to the QuantumToolbox.jl repository. It covers forking, cloning, setting up remotes, branching, committing, and pushing changes, as well as initiating a pull request.

```shell
git remote add qutip https://github.com/qutip/<repo-name>
git checkout main
git pull qutip main
git checkout -b <branch-name>
git add <file1> <file2> ...
git commit
git push -u origin <branch-name>
```

--------------------------------

### Initialize Makie plotting environment in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/extensions/makie.md

Demonstrates the necessary imports to enable Makie plotting functionality within QuantumToolbox.jl. Users must import both the main package and a specific Makie backend like CairoMakie.

```julia
using QuantumToolbox
using CairoMakie
```

--------------------------------

### Visualize Jaynes-Cummings Hamiltonian 3D Histogram in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Generates a 3D histogram visualization of the Jaynes-Cummings Hamiltonian using `matrix_histogram` from QuantumToolbox.jl. This provides a different perspective on the operator's structure. It requires `GLMakie` for optimal 3D rendering and `matrix_histogram`.

```julia
using QuantumToolbox
using QuantumOptics
using GLMakie # requires OpenGL: a backend better for 3D shading

# Assuming H is defined as in the first example
N = 5
a = tensor(destroy(N), qeye(2))
σm = tensor(qeye(N), destroy(2))
sx = tensor(qeye(N), sigmax())
H = a' * a + sx - 0.5 * (a * σm' + a' * σm)

fig = Figure(size = (500, 350))

matrix_histogram(H; location = fig[1,1])

fig
```

--------------------------------

### Converting Sparse Arrays to GPU Arrays

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/extensions/cuda.md

Illustrates the conversion of sparse CPU arrays (vectors and matrices) within QuantumObjects to GPU sparse arrays using the `cu` function. This includes examples for both `CuSparseVector` and `CuSparseMatrixCSC`, with an option to set the word size.

```julia
V = fock(2, 0; sparse=true) # CPU sparse vector
cu(V)
cu(V; word_size = 32)

M = sigmax() # CPU sparse matrix
cu(M)
cu(M; word_size = 32)
```

--------------------------------

### Make Commands for QuantumToolbox.jl Development

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/resources/contributing.md

This snippet lists common `make` commands used in the QuantumToolbox.jl project for development and maintenance. These commands simplify tasks such as setting up the environment, running tests, formatting code, building documentation, and generating changelogs.

```shell
make setup
make test
make format
make docs
make changelog
make help
```

--------------------------------

### Plot Wigner Function for Quantum States

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Calculates and plots the Wigner quasi-probability distribution for quantum states using QuantumToolbox.jl's plot_wigner function and Makie.jl. It visualizes the Wigner functions for coherent, thermal, and Fock states across a defined phase-space.

```julia
fig = Figure(size = (600, 200))

xvec = LinRange(-5, 5, 200)

_, ax1, _ = plot_wigner(ρ_coherent; xvec=xvec, yvec=xvec, location=fig[1,1])
ax1.title = "Coherent state"
ax1.aspect = 1 # make plot as square-size

_, ax2, _ = plot_wigner(ρ_thermal; xvec=xvec, yvec=xvec, location=fig[1,2])
ax2.title = "Thermal state"
ax2.aspect = 1 # make plot as square-size

_, ax3, _ = plot_wigner(ρ_fock; xvec=xvec, yvec=xvec, location=fig[1,3])
ax3.title = "Fock state"
ax3.aspect = 1 # make plot as square-size

fig
```

--------------------------------

### Accessing QobjEvo Properties

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/time_evolution/time_dependent.md

Demonstrates how to access various attributes and properties of a QuantumObjectEvolution (QobjEvo) object, such as its data, type, dimensions, and whether it represents a ket, bra, or operator.

```julia
using QuantumOptics

coef(p, t) = sin(π * t)
Ht = QobjEvo(sigmaz(), coef)

println(Ht.data)
println(Ht.type)
println(Ht.dims)
println(size(Ht))
println(shape(Ht)) # synonym of size(Ht)
println(length(Ht))
println(eltype(Ht)) # element type
println(isket(Ht)) # ket
println(isbra(Ht)) # bra
println(isoper(Ht)) # operator
println(isoperket(Ht)) # operator-ket
println(isoperbra(Ht)) # operator-bra
println(issuper(Ht)) # super operator
println(isconstant(Ht)) # time-independent or not
println(ishermitian(Ht)) # Hermitian
println(isherm(Ht)) # synonym of ishermitian(Ht)
println(issymmetric(Ht)) # symmetric
println(isposdef(Ht)) # positive definite (and Hermitian)
```

--------------------------------

### Visualize Jaynes-Cummings Hamiltonian Heatmap in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Generates a heatmap visualization of the Jaynes-Cummings Hamiltonian using QuantumToolbox.jl. This function takes the Hamiltonian operator as input and plots its matrix representation. It requires the `tensor`, `destroy`, `qeye`, `sigmax`, and `matrix_heatmap` functions from QuantumToolbox.jl and `Figure` from `Makie.jl`.

```julia
using QuantumToolbox
using QuantumOptics
using GLMakie

N = 5
a = tensor(destroy(N), qeye(2))
σm = tensor(qeye(N), destroy(2))
sx = tensor(qeye(N), sigmax())
H = a' * a + sx - 0.5 * (a * σm' + a' * σm)

fig = Figure(size = (500, 350))
matrix_heatmap(H; location = fig[1,1])

fig
```

--------------------------------

### Visualize Real and Imaginary Parts of Steady State Density Matrix in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/visualization.md

Visualizes the real and imaginary parts of a steady-state density matrix using `matrix_heatmap` from QuantumToolbox.jl. This allows for detailed inspection of the density matrix elements. It requires `sqrt`, `σm`, `a`, `steadystate`, `matrix_heatmap`, and `Figure`.

```julia
using QuantumToolbox
using QuantumOptics
using GLMakie

# Assuming H, a, σm are defined as in the previous example
N = 5
a = tensor(destroy(N), qeye(2))
σm = tensor(qeye(N), destroy(2))
sx = tensor(qeye(N), sigmax())
H = a' * a + sx - 0.5 * (a * σm' + a' * σm)

# collapse operators
c_ops = (
    sqrt(0.01) * σm,
    sqrt(0.01) * a,
)
ρss = steadystate(H, c_ops)

fig = Figure(size = (350, 500))

matrix_heatmap(ρss; method = Val(:real), location = fig[1,1])
matrix_heatmap(ρss; method = Val(:imag), location = fig[2,1])

fig
```

--------------------------------

### QuantumToolbox.jl: Basic steadystate Function Usage

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/steadystate.md

Demonstrates the fundamental ways to call the `steadystate` function in QuantumToolbox.jl. It accepts either a Hamiltonian and collapse operators, or a pre-defined Liouvillian. The output is the steady-state density matrix `ρ_ss`.

```julia
ρ_ss = steadystate(H)        # Hamiltonian
ρ_ss = steadystate(H, c_ops) # Hamiltonian and collapse operators
ρ_ss = steadystate(L)        # Liouvillian
```

--------------------------------

### Execute Quantum Simulations on GPU

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/README.md

Shows how to offload quantum dynamics calculations to a GPU by converting operators and states using the cu() function.

```julia
using QuantumToolbox
using CUDA
CUDA.allowscalar(false)

a_gpu = cu(destroy(N))
H_gpu = ω * a_gpu' * a_gpu
ψ0_gpu = cu(fock(N, 3))

c_ops = [sqrt(γ) * a_gpu]
e_ops = [a_gpu' * a_gpu]

sol = mesolve(H_gpu, ψ0_gpu, tlist, c_ops, e_ops = e_ops)
```

--------------------------------

### Calculate First-Order Optical Coherence Function g^(1)(τ) in Julia

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/two_time_corr_func.md

Calculates the first-order optical coherence function g^(1)(τ) for a system interacting with a thermal environment, starting from a coherent state. This function quantifies the coherence of the state over time. Requires `QuantumOptics.jl` and `QuantumToolbox.jl`.

```julia
using QuantumOptics, QuantumToolbox

τlist = LinRange(0, 10, 200)

# Hamiltonian
N = 15
a = destroy(N)
H = 2 * π * a' * a

# collapse operator
G1 = 0.75
n_th = 2.00  # bath temperature in terms of excitation number
c_ops = [
    sqrt(G1 * (1 + n_th)) * a,
    sqrt(G1 *      n_th)  * a'
]

# Initial state: coherent state
α = 2.5
ρ0 = coherent_dm(N, α)

# Calculate g^(1)(τ)
g1_tau = correlation_2op_1t(H, ρ0, τlist, c_ops, a', a; progress_bar = Val(false))

# The formula for g^(1)(τ) requires normalization using expectation values of a'a
# For a coherent state, <a'a> is constant, so we can use the initial value.
# For a general state, one would need to calculate <a'(t)a(t)> separately.
initial_na = expect(a' * a, ρ0)
normalized_g1_tau = g1_tau ./ sqrt.(initial_na .* initial_na)

# Plotting would typically follow using a plotting library like CairoMakie.jl
# For example:
# using CairoMakie
# fig = Figure()
# ax = Axis(fig[1, 1], xlabel = "τ", ylabel = "g^(1)(τ)")
# lines!(ax, τlist, abs.(normalized_g1_tau))
# fig
```

--------------------------------

### Access and Modify QuantumToolbox Global Settings

Source: https://github.com/qutip/quantumtoolbox.jl/blob/main/docs/src/users_guide/settings.md

Demonstrates how to inspect the current global configuration object and update specific parameters like tidyup_tol and auto_tidyup. These settings affect the behavior of various solver and utility functions within the package.

```julia
using QuantumToolbox

# Check current settings
QuantumToolbox.settings

# Modify settings
QuantumToolbox.settings.tidyup_tol = 1e-10
QuantumToolbox.settings.auto_tidyup = false

# Verify changes
QuantumToolbox.settings
```