diff --git a/src/QuantumOptics.jl b/src/QuantumOptics.jl
index 09ec19e4..6d166ab5 100644
--- a/src/QuantumOptics.jl
+++ b/src/QuantumOptics.jl
@@ -35,7 +35,9 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
         steadystate,
         timecorrelations,
         semiclassical,
-        stochastic
+        stochastic,
+        PauliBasis, PauliTransferMatrix, DensePauliTransferMatrix,
+        ChiMatrix, DenseChiMatrix
 
 
 include("sortedindices.jl")
@@ -86,6 +88,7 @@ module stochastic
     using .stochastic_definitions
 end
 include("printing.jl")
+include("pauli.jl")
 
 using .bases
 using .states
@@ -109,6 +112,7 @@ using .metrics
 using .spectralanalysis
 using .timecorrelations
 using .printing
+using .pauli
 
 
 end # module
diff --git a/src/pauli.jl b/src/pauli.jl
new file mode 100644
index 00000000..253cd3b3
--- /dev/null
+++ b/src/pauli.jl
@@ -0,0 +1,312 @@
+module pauli
+
+export PauliBasis, PauliTransferMatrix, DensePauliTransferMatrix,
+        ChiMatrix, DenseChiMatrix
+
+import Base: ==
+import Base: isapprox
+import Base: *
+
+using ..bases, ..spin, ..superoperators
+using ..operators: identityoperator, AbstractOperator
+using ..superoperators: SuperOperator
+using ..operators_dense: DenseOperator
+using ..spin: sigmax, sigmay, sigmaz
+using SparseArrays: sparse
+using LinearAlgebra: tr
+
+"""
+    PauliBasis(num_qubits::Int)
+
+Basis for an N-qubit space where `num_qubits` specifies the number of qubits.
+The dimension of the basis is 2²ᴺ.
+"""
+mutable struct PauliBasis{B<:Tuple{Vararg{Basis}}} <: Basis
+    shape::Vector{Int}
+    bases::B
+    function PauliBasis(num_qubits::Int)
+        shape = [2 for _ in 1:num_qubits]
+        bases = Tuple(SpinBasis(1//2) for _ in 1:num_qubits)
+        return new{typeof(bases)}(shape, bases)
+    end
+end
+==(pb1::PauliBasis, pb2::PauliBasis) = length(pb1.bases) == length(pb2.bases)
+
+"""
+    Base class for Pauli transfer matrix classes.
+"""
+abstract type PauliTransferMatrix{B1<:Tuple{PauliBasis, PauliBasis}, B2<:Tuple{PauliBasis, PauliBasis}} end
+
+
+"""
+    DensePauliTransferMatrix(B1, B2, data)
+
+DensePauliTransferMatrix stored as a dense matrix.
+"""
+mutable struct DensePauliTransferMatrix{B1<:Tuple{PauliBasis, PauliBasis},
+                                        B2<:Tuple{PauliBasis, PauliBasis},
+                                        T<:Matrix{Float64}} <: PauliTransferMatrix{B1, B2}
+    basis_l::B1
+    basis_r::B2
+    data::T
+    function DensePauliTransferMatrix(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{PauliBasis, PauliBasis},
+                                                                                BR<:Tuple{PauliBasis, PauliBasis},
+                                                                                T<:Matrix{Float64}}
+        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
+           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
+            throw(DimensionMismatch())
+        end
+        new{BL, BR, T}(basis_l, basis_r, data)
+    end
+end
+
+PauliTransferMatrix(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis} = ptm
+
+function *(ptm0::DensePauliTransferMatrix{B, B, Matrix{Float64}},
+           ptm1::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    return DensePauliTransferMatrix(ptm0.basis_l, ptm1.basis_r, ptm0.data*ptm1.data)
+end
+
+"""
+    Base class for χ (process) matrix classes.
+"""
+abstract type ChiMatrix{B1<:Tuple{PauliBasis, PauliBasis}, B2<:Tuple{PauliBasis, PauliBasis}} end
+
+"""
+    DenseChiMatrix(b, b, data)
+
+DenseChiMatrix stored as a dense matrix.
+"""
+mutable struct DenseChiMatrix{B1<:Tuple{PauliBasis, PauliBasis},
+                              B2<:Tuple{PauliBasis, PauliBasis},
+                              T<:Matrix{ComplexF64}} <: PauliTransferMatrix{B1, B2}
+    basis_l::B1
+    basis_r::B2
+    data::T
+    function DenseChiMatrix(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{PauliBasis, PauliBasis},
+                                                                                BR<:Tuple{PauliBasis, PauliBasis},
+                                                                                T<:Matrix{ComplexF64}}
+        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
+           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
+            throw(DimensionMismatch())
+        end
+        new{BL, BR, T}(basis_l, basis_r, data)
+    end
+end
+
+ChiMatrix(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis} = chi_matrix
+
+"""
+A dictionary that represents the Pauli algebra - for a pair of Pauli operators
+σᵢσⱼ information about their product is given under the key "ij". The first
+element of the dictionary value is the Pauli operator, and the second is the
+scalar multiplier. For example, σ₀σ₁ = σ₁, and `"01" => ("1", 1)`.
+"""
+const pauli_multiplication_dict = Dict(
+  "00" => ("0", 1.0+0.0im),
+  "23" => ("1", 0.0+1.0im),
+  "30" => ("3", 1.0+0.0im),
+  "22" => ("0", 1.0+0.0im),
+  "21" => ("3", -0.0-1.0im),
+  "10" => ("1", 1.0+0.0im),
+  "31" => ("2", 0.0+1.0im),
+  "20" => ("2", 1.0+0.0im),
+  "01" => ("1", 1.0+0.0im),
+  "33" => ("0", 1.0+0.0im),
+  "13" => ("2", -0.0-1.0im),
+  "32" => ("1", -0.0-1.0im),
+  "11" => ("0", 1.0+0.0im),
+  "03" => ("3", 1.0+0.0im),
+  "12" => ("3", 0.0+1.0im),
+  "02" => ("2", 1.0+0.0im),
+)
+
+"""
+    multiply_pauli_matirices(i4::String, j4::String)
+
+A function to algebraically determine result of multiplying two
+(N-qubit) Pauli matrices. Each Pauli matrix is represented by a string
+in base 4. For example, σ₃⊗σ₀⊗σ₂ would be "302". The product of any pair of
+Pauli matrices will itself be a Pauli matrix multiplied by any of the 1/4 roots
+of 1.
+"""
+cache_multiply_pauli_matrices() = begin
+    local pauli_multiplication_cache = Dict()
+    function _multiply_pauli_matirices(i4::String, j4::String)
+        if (i4, j4) ∉ keys(pauli_multiplication_cache)
+            pauli_multiplication_cache[(i4, j4)] = mapreduce(x -> pauli_multiplication_dict[prod(x)],
+                                                             (x,y) -> (x[1] * y[1], x[2] * y[2]),
+                                                             zip(i4, j4))
+        end
+        return pauli_multiplication_cache[(i4, j4)]
+    end
+end
+multiply_pauli_matirices = cache_multiply_pauli_matrices()
+
+function *(chi_matrix0::DenseChiMatrix{B, B, Matrix{ComplexF64}},
+           chi_matrix1::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+
+    num_qubits = length(chi_matrix0.basis_l[1].shape)
+    sop_dim = 2 ^ prod(chi_matrix0.basis_l[1].shape)
+    ret = zeros(ComplexF64, (sop_dim, sop_dim))
+
+    for ijkl in Iterators.product(0:(sop_dim-1),
+                                  0:(sop_dim-1),
+                                  0:(sop_dim-1),
+                                  0:(sop_dim-1))
+        i, j, k, l = ijkl
+        if (chi_matrix0.data[i+1, j+1] != 0.0) & (chi_matrix1.data[k+1, l+1] != 0.0)
+            i4, j4, k4, l4 = map(x -> string(x, base=4, pad=2), ijkl)
+
+            pauli_product_ik = multiply_pauli_matirices(i4, k4)
+            pauli_product_lj = multiply_pauli_matirices(l4, j4)
+
+            ret[parse(Int, pauli_product_ik[1], base=4)+1,
+                parse(Int, pauli_product_lj[1], base=4)+1] += (pauli_product_ik[2] * pauli_product_lj[2] * chi_matrix0.data[i+1, j+1] * chi_matrix1.data[k+1, l+1])
+        end
+    end
+    return DenseChiMatrix(chi_matrix0.basis_l, chi_matrix0.basis_r, ret / 2^num_qubits)
+end
+
+
+# TODO MAKE A GENERATOR FUNCTION
+"""
+    pauli_operators(num_qubits::Int)
+
+Generate a list of N-qubit Pauli operators.
+"""
+function pauli_operators(num_qubits::Int)
+    pauli_funcs = (identityoperator, sigmax, sigmay, sigmaz)
+    po = []
+    for paulis in Iterators.product((pauli_funcs for _ in 1:num_qubits)...)
+        basis_vector = reduce(⊗, f(SpinBasis(1//2)) for f in paulis)
+        push!(po, basis_vector)
+    end
+    return po
+end
+
+"""
+    pauli_basis_vectors(num_qubits::Int)
+
+Generate a matrix of basis vectors in the Pauli representation given a number
+of qubits.
+"""
+function pauli_basis_vectors(num_qubits::Int)
+    po = pauli_operators(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    return mapreduce(x -> sparse(reshape(x.data, sop_dim)), (x, y) -> [x y], po)
+end
+
+"""
+    PauliTransferMatrix(sop::DenseSuperOperator)
+
+Convert a superoperator to its representation as a Pauli transfer matrix.
+"""
+function PauliTransferMatrix(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(sop.basis_l[1].bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    data = Matrix{Float64}(undef, (sop_dim, sop_dim))
+    data .= real.(pbv' * sop.data * pbv / √sop_dim)
+    return DensePauliTransferMatrix(sop.basis_l, sop.basis_r, data)
+end
+
+SuperOperator(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis = spre(unitary) * spost(unitary')
+SuperOperator(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis} = sop
+
+"""
+    SuperOperator(ptm::DensePauliTransferMatrix)
+
+Convert a Pauli transfer matrix to its representation as a superoperator.
+"""
+function SuperOperator(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(ptm.basis_l[1].bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    sop_dim = 4 ^ num_qubits
+    data = Matrix{ComplexF64}(undef, (sop_dim, sop_dim))
+    data .= pbv * ptm.data * pbv' / √sop_dim
+    return DenseSuperOperator(ptm.basis_l, ptm.basis_r, data)
+end
+
+"""
+    PauliTransferMatrix(unitary::DenseOperator)
+
+Convert an operator, presumably a unitary operator, to its representation as a
+Pauli transfer matrix.
+"""
+PauliTransferMatrix(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis = PauliTransferMatrix(SuperOperator(unitary))
+
+"""
+    ChiMatrix(unitary::DenseOperator)
+
+Convert an operator, presumably a unitary operator, to its representation as a χ matrix.
+"""
+function ChiMatrix(unitary::DenseOperator{B, B, Matrix{ComplexF64}}) where B <: PauliBasis
+    num_qubits = length(unitary.basis_l.bases)
+    pbv = pauli_basis_vectors(num_qubits)
+    aj = pbv' * reshape(unitary.data, 4 ^ num_qubits)
+    return DenseChiMatrix((unitary.basis_l, unitary.basis_l), (unitary.basis_r, unitary.basis_r), aj * aj' / (2 ^ num_qubits))
+end
+
+"""
+    ChiMatrix(sop::DenseSuperOperator)
+
+Convert a superoperator to its representation as a Chi matrix.
+"""
+function ChiMatrix(sop::DenseSuperOperator{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(sop.basis_l)
+    sop_dim = 4 ^ num_qubits
+    po = pauli_operators(num_qubits)
+    data = Matrix{ComplexF64}(undef, (sop_dim, sop_dim))
+    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
+        data[idx, jdx] = tr((spre(po[idx]) * spost(po[jdx])).data' * sop.data) / √sop_dim
+    end
+    return DenseChiMatrix(sop.basis_l, sop.basis_r, data)
+end
+
+"""
+    PauliTransferMatrix(chi_matrix::DenseChiMatrix)
+
+Convert a χ matrix to its representation as a Pauli transfer matrix.
+"""
+function PauliTransferMatrix(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    num_qubits = length(chi_matrix.basis_l)
+    sop_dim = 4 ^ num_qubits
+    po = pauli_operators(num_qubits)
+    data = Matrix{Float64}(undef, (sop_dim, sop_dim))
+    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
+        data[idx, jdx] = tr(mapreduce(x -> po[idx] * po[x[1]] * po[jdx] * po[x[2]] * chi_matrix.data[x[1], x[2]],
+                                      +,
+                                      Iterators.product(1:16, 1:16)).data) / sop_dim |> real
+    end
+    return DensePauliTransferMatrix(chi_matrix.basis_l, chi_matrix.basis_r, data)
+end
+
+"""
+    SuperOperator(chi_matrix::DenseChiMatrix)
+
+Convert a χ matrix to its representation as a superoperator.
+"""
+function SuperOperator(chi_matrix::DenseChiMatrix{B, B, Matrix{ComplexF64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    return SuperOperator(PauliTransferMatrix(chi_matrix))
+end
+
+"""
+    ChiMatrix(ptm::DensePauliTransferMatrix)
+
+Convert a Pauli transfer matrix to its representation as a χ matrix.
+"""
+function ChiMatrix(ptm::DensePauliTransferMatrix{B, B, Matrix{Float64}}) where B <: Tuple{PauliBasis, PauliBasis}
+    return ChiMatrix(SuperOperator(ptm))
+end
+
+"""Equality for all varieties of superoperators."""
+==(sop1::T, sop2::T) where T<:Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix} = sop1.data == sop2.data
+==(sop1::Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix}, sop2::Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix}) = false
+
+"""Approximate equality for all varieties of superoperators."""
+function isapprox(sop1::T, sop2::T; kwargs...) where T<:Union{DensePauliTransferMatrix, DenseSuperOperator, DenseChiMatrix}
+    return isapprox(sop1.data, sop2.data; kwargs...)
+end
+
+end # end module
diff --git a/test/runtests.jl b/test/runtests.jl
index 08db82cd..ce75d3f9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,6 +46,7 @@ names = [
     "test_stochastic_schroedinger.jl",
     "test_stochastic_master.jl",
     "test_stochastic_semiclassical.jl",
+    "test_pauli.jl",
 
     "test_printing.jl"
 ]
diff --git a/test/test_pauli.jl b/test/test_pauli.jl
new file mode 100644
index 00000000..66c51944
--- /dev/null
+++ b/test/test_pauli.jl
@@ -0,0 +1,92 @@
+using LinearAlgebra
+using Test
+
+using QuantumOptics
+
+@testset "pauli" begin
+
+@test_throws MethodError PauliBasis(1.4)
+
+# Test conversion of unitary matrices to superoperators.
+q2 = PauliBasis(2)
+q3 = PauliBasis(3)
+CZ = DenseOperator(q2, q2, diagm(0 => [1,1,1,-1]))
+CZ_sop = SuperOperator(CZ)
+
+# Test conversion of unitary matrices to superoperators.
+@test diag(CZ_sop.data) ==  ComplexF64[1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,-1,1]
+@test CZ_sop.basis_l == CZ_sop.basis_r == (q2, q2)
+
+# Test conversion of superoperator to Pauli transfer matrix.
+CZ_ptm = PauliTransferMatrix(CZ_sop)
+
+# Test DensePauliTransferMatrix constructor.
+@test_throws DimensionMismatch DensePauliTransferMatrix((q2, q2), (q3, q3), CZ_ptm.data)
+@test DensePauliTransferMatrix((q2, q2), (q2, q2), CZ_ptm.data) == CZ_ptm
+
+@test all(isapprox.(CZ_ptm.data[[1,30,47,52,72,91,117,140,166,185,205,210,227,256]], 1))
+@test all(isapprox.(CZ_ptm.data[[106,151]], -1))
+
+@test CZ_ptm == PauliTransferMatrix(ChiMatrix(CZ))
+
+# Test construction of non-symmetric unitary.
+CNOT = DenseOperator(q2, q2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
+CNOT_sop = SuperOperator(CNOT)
+CNOT_chi = ChiMatrix(CNOT)
+CNOT_ptm = PauliTransferMatrix(CNOT)
+
+@test CNOT_sop.basis_l == CNOT_sop.basis_r == (q2, q2)
+@test CNOT_chi.basis_l == CNOT_chi.basis_r == (q2, q2)
+@test CNOT_ptm.basis_l == CNOT_ptm.basis_r == (q2, q2)
+
+@test all(isapprox.(imag.(CNOT_sop.data), 0))
+@test all(isapprox.(imag.(CNOT_chi.data), 0))
+@test all(isapprox.(imag.(CNOT_ptm.data), 0))
+
+@test all(isapprox.(CNOT_sop.data[[1,18,36,51,69,86,104,119,141,158,176,191,201,218,236,251]], 1))
+@test all(isapprox.(CNOT_chi.data[[1,2,13,17,18,29,193,194,205,222]], 1))
+@test all(isapprox.(CNOT_chi.data[[14,30,206,209,210,221]], -1))
+@test all(isapprox.(CNOT_ptm.data[[1,18,47,64,70,85,108,138,153,183,205,222,227,244,]], 1))
+@test all(isapprox.(CNOT_ptm.data[[123,168]], -1))
+
+# Test DenseChiMatrix constructor.
+@test_throws DimensionMismatch DenseChiMatrix((q2, q2), (q3, q3), CNOT_chi.data)
+@test DenseChiMatrix((q2, q2), (q2, q2), CNOT_chi.data) == CNOT_chi
+
+# Test equality and conversion among all three bases.
+ident = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
+
+IDENT = DenseOperator(q2, ident)
+
+IDENT_sop = SuperOperator(IDENT)
+IDENT_chi = ChiMatrix(IDENT)
+IDENT_ptm = PauliTransferMatrix(IDENT)
+
+@test ChiMatrix(IDENT_sop) == IDENT_chi
+@test ChiMatrix(IDENT_ptm) == IDENT_chi
+@test SuperOperator(IDENT_chi) == IDENT_sop
+@test SuperOperator(IDENT_ptm) == IDENT_sop
+@test PauliTransferMatrix(IDENT_sop) == IDENT_ptm
+@test PauliTransferMatrix(IDENT_chi) == IDENT_ptm
+
+# Test approximate equality and conversion among all three bases.
+cphase = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 exp(1im*.6)]
+
+CPHASE = DenseOperator(q2, cphase)
+
+CPHASE_sop = SuperOperator(CPHASE)
+CPHASE_chi = ChiMatrix(CPHASE)
+CPHASE_ptm = PauliTransferMatrix(CPHASE)
+
+@test isapprox(ChiMatrix(CPHASE_sop), CPHASE_chi)
+@test isapprox(ChiMatrix(CPHASE_ptm), CPHASE_chi)
+@test isapprox(SuperOperator(CPHASE_chi), CPHASE_sop)
+@test isapprox(SuperOperator(CPHASE_ptm), CPHASE_sop)
+@test isapprox(PauliTransferMatrix(CPHASE_sop), CPHASE_ptm)
+@test isapprox(PauliTransferMatrix(CPHASE_chi), CPHASE_ptm)
+
+# Test composition.
+@test isapprox(ChiMatrix(CPHASE) * ChiMatrix(CNOT), ChiMatrix(CPHASE * CNOT))
+@test isapprox(PauliTransferMatrix(CPHASE) * PauliTransferMatrix(CNOT), PauliTransferMatrix(CPHASE * CNOT))
+
+end # testset