JuliaでPauli行列のコンストラクタをつくる2

以前の私の記事では、Pauli行列を"行列のまま"保持してリスト化したコンストラクタを作成しました。

今回はそれらを発展させて、Pauli行列自体もコンストラクタにしようと試みました。

コード

次のようなコードを作成しました。

const BASE_PAULI = [
    Matrix{ComplexF64}([1. 0.; 0. 1.]),
    Matrix{ComplexF64}([0. 1.; 1. 0.]),
    Matrix{ComplexF64}([0. -im; im 0.]),
    Matrix{ComplexF64}([1. 0.; 0. -1.]),
]

abstract type AbstractPauli{N} <: AbstractArray{ComplexF64,N} end

struct PauliMatrix <: AbstractPauli{2}
  matrix::Matrix
  space::String
end
Base.:+(M::PauliMatrix,N::PauliMatrix) = (
  if M.space == N.space
    PauliMatrix(M.matrix + N.matrix, M.space)
  else
    error("Pauli matrices M ⊂ $(M.space) and N ⊂ $(N.space) mismatch.")
  end
  )
Base.:-(M::PauliMatrix) = PauliMatrix(-M.matrix,M.space)
Base.:-(M::PauliMatrix,N::PauliMatrix) = (
  if M.space == N.space
    PauliMatrix(M.matrix - N.matrix, M.space)
  else
    error("Pauli matrices M ⊂ $(M.space) and N ⊂ $(N.space) mismatch.")
  end
  )
Base.:*(M::PauliMatrix,N::PauliMatrix) = (
  if M.space == N.space
    PauliMatrix(M.matrix * N.matrix, M.space)
  else
    PauliMatrix(Base.kron(M.matrix,N.matrix),"$(M.space)⊗$(N.space)")
  end
  )
Base.:*(a::Number,M::PauliMatrix) = PauliMatrix( a * M.matrix, M.space)
Base.:*(M::PauliMatrix,a::Number) = PauliMatrix( a * M.matrix, M.space)
Base.kron(M::PauliMatrix,N::PauliMatrix) = (
    PauliMatrix(Base.kron(M.matrix,N.matrix),"$(M.space)⊗$(N.space)")
  )
⊗(M::PauliMatrix,N::PauliMatrix) = Base.kron(M,N)


Base.@kwdef struct PauliMatrices <:AbstractPauli{1}
    space::String = "_"
    
    I = PauliMatrix(BASE_PAULI[1],space)
    x = PauliMatrix(BASE_PAULI[2],space)
    y = PauliMatrix(BASE_PAULI[3],space)
    z = PauliMatrix(BASE_PAULI[4],space)

    matrix = Vector([I,x,y,z])
end
PauliMatrices(space) = PauliMatrices(space=space)

Base.getindex(P::PauliMatrices, ind::Integer) = P.matrix[ind+1]
Base.getindex(P::PauliMatrix, i::Integer, j::Integer) = P.matrix[i,j]
Base.size(P::AbstractPauli) = Base.size(P.matrix)
Base.show(P::AbstractPauli) = Base.show(P.matrix)
Base.print_array(io::IO, P::AbstractPauli) = Base.print_matrix(io, P.matrix)
Base.array_summary(
  io::IO,
  P::AbstractPauli,
  inds::Tuple{Vararg{Base.OneTo}}
  ) = begin
    print(io, Base.dims2string(length.(inds)), " ")
    Base.showarg(io, P, true)
    if length(P.space) < min(18,displaysize(io)[2])
      print(io, "($(P.space))")
    else
      space_list = split(P.space,"⊗")
      print(io, string(P.space)[1:10],"⋅⋅⋅⊗$(space_list[end]))")
    end
end

具体例

簡単に説明していきます。

• Pauli行列の作成
  Pauli行列(のカタマリ)は次のように作成します。

  σ = PauliMatrices("spin")
  τ = PauliMatrices("orbital")
  

  ここで、spaceと呼んでいるのは、Pauli行列の表している空間です。例えば今回で言えばσはスピンの、τは軌道の自由度を表しています。
  spaceを入力しない場合は、"_"という空間になるように初期化しています。

• Pauli行列へのアクセス
  こちらは以前と同様の方法でアクセスできます。

  julia> σ[0]
  2×2 PauliMatrix(spin):
   1.0+0.0im  0.0+0.0im
   0.0+0.0im  1.0+0.0im
  
  julia> σ[1]
  2×2 PauliMatrix(spin):
   0.0+0.0im  1.0+0.0im
   1.0+0.0im  0.0+0.0im
  
  julia> σ.y
  2×2 PauliMatrix(spin):
   0.0+0.0im  0.0-1.0im
   0.0+1.0im  0.0+0.0im
  

  ただし、以前と違う点として、属性がPauliMatrixという新たに定義したものに変わっています

• Pauli行列の演算
  今回わざわざPauli行列をコンストラクタで定義した理由がここです。
  異なる空間でPauli行列を積演算する場合には、クロネッカー積で2つのPauli行列を計算する必要があります。
  $$
  \tau_{\mu}\sigma_{\nu} = \tau_{\mu}\otimes\sigma_{\nu}
  $$
  これを自動で判別してくれるように変更しました。
  以下具体例です。

  julia> τ[1]*τ[2]
  2×2 PauliMatrix(orbital):
   0.0+1.0im  0.0+0.0im
   0.0+0.0im  0.0-1.0im
  
  julia> τ[1]*σ[2]
  4×4 PauliMatrix(orbital⊗spin):
   0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0-1.0im
   0.0+0.0im  0.0+0.0im  0.0+1.0im  0.0+0.0im
   0.0+0.0im  0.0-1.0im  0.0+0.0im  0.0+0.0im
   0.0+1.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  

  このように、同じ*という積を用いても、自動で空間を判別して計算手法を判別できるように変更しました。

  また、spaceも自動で次のように変更されます。

  julia> (τ[1]*σ[2]).space
  "orbital⊗spin"
  

  逆に、空間を定義していない（"undefined"空間の）行列では、インスタンス名が異なっていても、同等として扱うので注意してください。

  julia> σ_ = PauliMatrices();
  julia> τ_ = PauliMatrices();
  
  julia> τ_[1]*σ_[2]
  2×2 PauliMatrix(_):
   0.0+1.0im  0.0+0.0im
   0.0+0.0im  0.0-1.0im
  

  一方、⊗を用いれば、同空間においてもクロネッカー積が施せます。

  julia> τ_[1]⊗σ_[2]
  4×4 PauliMatrix(_⊗_):
   0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0-1.0im
   0.0+0.0im  0.0+0.0im  0.0+1.0im  0.0+0.0im
   0.0+0.0im  0.0-1.0im  0.0+0.0im  0.0+0.0im
   0.0+1.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  

  和演算については、同空間でないときにエラーを吐くようにしました。

  julia> τ[1] + σ[2]
  ERROR: The definition spaces of the respective Pauli matrices mismatch
  Stacktrace:
   [1] error(s::String)
     @ Base .\error.jl:35
   [2] +(M::PauliMatrix, N::PauliMatrix)
     @ Main R:\\pauli.jl:13
   [3] top-level scope
     @ REPL[18]:1
   [4] top-level scope
  

さいごに

まだまだ必要になるPauli行列の特性は多くありそうなので（expの演算とかは定義する価値があったりするかもしれないです）、いろいろと実装できたらと思います。
皆様からのご提案などもいただけたら非常にうれしいです。