Public release

Q3/Documentation/English/SearchIndex/5/indexMetadata.wl

@@ -1,6 +1,6 @@
 <|"IndexedPaths" -> None, "Driver" -> "Lucene", 
- "CreationDate" -> DateObject[{2024, 6, 17, 21, 22, 
-    26.261208`8.171889681570036}, "Instant", "Gregorian", 9.], 
+ "CreationDate" -> DateObject[{2024, 6, 18, 19, 29, 
+    29.949547`8.228965239131098}, "Instant", "Gregorian", 9.], 
  "Version" -> 5, "Synonyms" -> None, "Method" -> "BM25", 
  "Language" -> "English", "ContentFieldOptions" -> 
   <|"Title" -> <|"Stored" -> True, "Weight" -> 2|>, 

Q3/Kernel/Abby.wl

@@ -314,6 +314,8 @@ PseudoDivide[x_, 0. + 0. I] = 0. + 0. I
 PseudoDivide[x_, y_] := x/y
 
 
+(**** <ZeroQ> ****)
+
 ZeroQ::usage = "ZeroQ[x] returns True if x approximately equals to zero.\nZeroQ[x, \[Delta]] returns True if |x| \[LessEqual] \[Delta]."
 
 ZeroQ::tolnn = "Tolerence specification `` must be a non-negative number."
@@ -325,7 +327,9 @@ ZeroQ[x_, del_?NonNegative] := TrueQ[Chop[x, del] == 0]
 ZeroQ[x_, del_] := (
   Message[ZeroQ::tolnn, del];
   ZeroQ[x, Abs @ del]
- )
+)
+
+(**** </ZeroQ> ****)
 
 
 (**** <CountsFor> ****)

Q3/Kernel/Abel.wl

@@ -84,7 +84,6 @@ $::usage = "$ is a flavor index referring to the species itself."
   CliffordGroupElements, FullCliffordGroupElements,
   CliffordQ };
 
-
 (**** <compatibiliy> ****)
 
 (* Now comes as a built-in function with v13.1, but with an additional

Q3/Kernel/Gottesman.wl

@@ -1,4 +1,3 @@
-0(* -*- mode:math -*- *)
 (* stabilizer formalism *)
 BeginPackage["Q3`"]
 
@@ -13,7 +12,9 @@ BeginPackage["Q3`"]
 
 { CliffordDecompose = CliffordFactor };
 
-{ Stabilizer, StabilizerStateQ, StabilizerStateCount };
+{ Stabilizer, StabilizerGenerators };
+
+{ StabilizerStateQ, StabilizerStateCount };
 
 { GottesmanGroup = BinarySymplecticGroup,
   GottesmanGroupElements = BinarySymplecticGroupElements };
@@ -552,40 +553,31 @@ GottesmanMerge[xx_?MatrixQ, zz_?MatrixQ] := MapThread[Riffle, {xx, zz}]
 
 (**** <Stabilizer> ****)
 
-Stabilizer::usage = "Stabilzier[state] returns the stabilizer subgroup of the Pauli group that stabilizes state, which may be a column vecotr or expressed in terms of Ket[\[Ellipsis]] or Ket[<|\[Ellipsis]|>].\nStabilizer[state, {s1,s2,\[Ellipsis]}] assumes that state belongs to the Hilbert space associated with qubits {s1,s2,\[Ellipsis]}.\nStabilizer[graph] returns the stabilizer subgroup of the Pauli group that stabilizes the graph state associated with  graph.\nStabilizer[graph, vtx] gives the operator associated with the vertex vtx that stabilize the graph state associated with graph."
+Stabilizer::usage = "Stabilzier[state] returns the stabilizer subgroup of the Pauli group that stabilizes state, which may be a column vecotr or expressed in terms of Ket[\[Ellipsis]] or Ket[<|\[Ellipsis]|>].\nStabilizer[state, {s1,s2,\[Ellipsis]}] assumes that state belongs to the Hilbert space associated with qubits {s1,s2,\[Ellipsis]}.\nStabilizer[graph] returns the stabilizer subgroup of the Pauli group that stabilizes the graph state associated with graph.\nStabilizer[graph, vtx] returns the operator associated with vertex vtx (the so-called correlation operator on vtx) that stabilizes the graph state associated with graph."
 
 Stabilizer::notss = "`` is not a stabilizer state."
 
 Stabilizer[vec_?VectorQ] := With[
   { mm = getStabilizer[vec] },
-  If[ FailureQ[mm],
-    Message[Stabilizer::notss, expr];
-    Return[{}]
-   ];
+  If[FailureQ[mm], Message[Stabilizer::notss, expr]; Return @ {}];
   Map[ThePauli, Keys @ mm] * Values[mm]
- ] /; IntegerQ @ Log[2, Length @ vec]
+] /; IntegerQ @ Log[2, Length @ vec]
 
 
 Stabilizer[expr_] := With[
   { mm = getStabilizer @ Matrix[expr] },
-  If[ FailureQ[mm],
-    Message[Stabilizer::notss, expr];
-    Return[{}]
-   ];
+  If[FailureQ[mm], Message[Stabilizer::notss, expr]; Return @ {}];
   Map[Pauli, Keys @ mm] * Values[mm]
- ] /; fPauliKetQ[expr]
+] /; fPauliKetQ[expr]
 
 
 Stabilizer[expr_] := Stabilizer[expr, Qubits @ expr]
 
 Stabilizer[expr_, ss:{__?QubitQ}] := With[
   { mm = getStabilizer @ Matrix[expr, ss] },
-  If[ FailureQ[mm],
-    Message[Stabilizer::notss, expr];
-    Return[{}]
-   ];
+  If[FailureQ[mm], Message[Stabilizer::notss, expr]; Return @ {}];
   Elaborate[ (Multiply @@@ FlavorThread[ss, Keys @ mm]) * Values[mm] ]
- ]
+]
 
 
 getStabilizer[vec_?VectorQ] := Module[
@@ -594,7 +586,7 @@ getStabilizer[vec_?VectorQ] := Module[
   tsr = PauliDecompose[rho];
   If[Not[Equal @@ Abs[Values @ tsr]], Return @ $Failed];
   Sign /@ tsr
- ]
+]
 
 
 Stabilizer[g_Graph] := Stabilizer @ GraphState[g]
@@ -607,8 +599,31 @@ Stabilizer[g_Graph, vtx_?QubitQ] := Module[
     adj },
   adj = AdjacencyList[g, new|new[$]];
   vtx[1] ** Apply[Multiply, Through[adj[3]]]
- ]
+]
+
+(**** </Stabilizer> ****)
+
+
+(**** <StabilizerGenerators> ****)
+
+StabilizerGenerators::usage = "StabilizerGenerators[{g1, g2, \[Ellipsis]}] returns a list of independent generators of stabilizer subgroup {g1, g2, \[Ellipsis]} of the Pauli group."
 
+StabilizerGenerators[grp_List] := Module[
+  { ss = Qubits[grp],
+    gg, cc
+  },
+  gg = GottesmanVector[#, ss]& /@ grp;
+  gg = Orthogonalize[gg, Mod[Dot[#1, #2], 2]&];
+  gg = DeleteCases[gg, {0..}];
+  gg = FromGottesmanVector[#, ss]& /@ gg;
+  cc = FirstCase[Coefficient[grp, #], Except[0]] & /@ gg;
+  cc * gg
+]
+
+(**** </StabilizerGenerators> ****)
+
+
+(**** <StabilizerState> ****)
 
 StabilizerStateQ::usage = "StabilizerStateQ[state] returns True if state is a stabilizer state, a state that can be stabilized by a non-trivial subgroup of the Pauli group; and False otherwise. The state may be a column vector or expressed in terms of Ket[\[Ellipsis]] or Ket[<|\[Ellipsis]|>]."
 
@@ -617,7 +632,7 @@ StabilizerStateQ::notqbt = "`` is not a state vector for qubits."
 StabilizerStateQ[vec_?VectorQ] := (
   Message[StabilizerStateQ::notqbt, vec];
   False
- ) /; Not @ IntegerQ @ Log[2, Length @ vec]
+) /; Not @ IntegerQ @ Log[2, Length @ vec]
 
 StabilizerStateQ[vec_?VectorQ] :=
   TrueQ[Equal @@ Abs[Values @ PauliDecompose @ Dyad[vec, vec]]]
@@ -636,7 +651,7 @@ StabilizerStateCount[n_Integer] :=
 StabilizerStateCount[ss:{___?QubitQ}] :=
   StabilizerStateCount[Length @ ss]
 
-(**** </Stabilizer> ****)
+(**** </StabilizerState> ****)
 
 
 (**** <GottesmanGroup> ****)

Q3/Kernel/Wick.wl

@@ -7,57 +7,49 @@ BeginPackage["Q3`"]
 
 { WickRandomCircuit };
 
-Begin["`Private`"]
-
-(**** <WickLogarithmicNegtivity> ****)
-
-(* See also: Shapourian and Ryu (2017, 2019) *)
-
-WickLogarithmicNegativity::usage = "WickLogarithmicNegativity[{grn,anm}, {k1,k2,\[Ellipsis]}] returns the fermionic negativity over modes k1, k2, \[Ellipsis] of the fermionic Gaussian states characterized by the Green's function grn (in an L\[Times]L matrix for L fermion modes) and anomalous Green's function anm (also in an L\[Times]L matrix).\nWickLogarithmicNegativity[grn, {k1,k2,\[Ellipsis]}] is equivalent to WickLogarithmicNegativity[{grn,0}, {k1,k2,\[Ellipsis]}]."
-
-WickLogarithmicNegativity[kk:{__Integer}][any_] :=
-  WickLogarithmicNegativity[any, kk]
-
-WickLogarithmicNegativity[grn_?MatrixQ, kk:{__Integer}] :=
-  WickTimeReversalMoment[1/2, grn, kk]
-
-WickLogarithmicNegativity[{grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}] := 
-  WickTimeReversalMoment[1/2, {grn, anm}, kk]
-
-(**** </WickLogarithmicNegtivity> ****)
+{ EntanglementEntropy }; (* mainly in vonNeumann.wl *)
 
+Begin["`Private`"]
 
 (**** <WickTimeReversalMoment> ****)
 
 (* See also: Shapourian and Ryu (2017, 2019) *)
 
 WickTimeReversalMoment::usage = "WickTimeReversalMoment[\[Alpha], {grn,anm}, {k1,k2,\[Ellipsis]}] returns the \[Alpha]th moment of partial time reversal over the fermion modes (species) k1, k2, \[Ellipsis] for the fermionic Gaussian state characterized by the Green's funciton grn (in an L\[Times]L matrix for L fermion modes) and anomalous Green's function anm (also in an L\[Times]L matrix).\nWickTimeReversalMoment[\[Alpha], grn, {k1,k2,\[Ellipsis]}] is equivalent to WickTimeReversalMoment[\[Alpha], {grn,0}, {k1,k2,\[Ellipsis]}]."
 
-WickTimeReversalMoment::sing = "Matrix `` is singular; PseudoInverse is used instead of Inverse."
+WickTimeReversalMoment::sing = "The matrix is tamed according to option \"Epsilon\"."
+
+Options[WickTimeReversalMoment] = { "Epsilon" -> 1.25*^-5 } (* 1.0*^-8 is also a reasonable choice *)
 
 
-WickTimeReversalMoment[alpha_, grn_?MatrixQ, kk:{__Integer}] :=
-  WickTimeReversalMoment[alpha, {grn, Zero[Dimensions @ grn]}, kk]
+WickTimeReversalMoment[alpha_, grn_?MatrixQ, kk:{__Integer}, opts___?OptionQ] :=
+  WickTimeReversalMoment[alpha, {grn, Zero[Dimensions @ grn]}, kk, opts]
 
-WickTimeReversalMoment[alpha_, {grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}] := Module[
+WickTimeReversalMoment[alpha_, {grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}, OptionsPattern[]] := Module[
   { dd = Length[grn],
-    gg, zr, id, xx, zz, uu, ww, pf, dgn, off
+    gg, oo, id, xx, zz, uu, ww, pf, dgn, off
   },
-  zr = Zero[dd, dd];
+  oo = Zero[dd, dd];
   id = One[dd];
   xx = CircleTimes[ThePauli[1], id];
   zz = CircleTimes[ThePauli[3], id];
-  
+
   (* \Gamma *)
+  gg = id - grn;
+  Check[
+    Inverse @ gg,
+    Message[WickTimeReversalMoment::sing];
+    gg = tameMatrix[gg, OptionValue["Epsilon"]]
+  ];
   gg = ArrayFlatten @ {
-    {id - grn, anm},
-    {Topple @ anm, Transpose[grn] - id}
+    {gg, anm},
+    {Topple @ anm, -Transpose[gg]}
   };
-  pf = Check[Inverse[gg.xx], PseudoInverse[gg.xx]];
-  pf = Power[Pfaffian[pf], 2];
+
+  pf = Power[Pfaffian[Inverse[gg.xx]], 2];
 
   (* \Omega of the density operator \rho *)
-  ww = Check[Inverse[gg] - zz, PseudoInverse[gg] - zz];
+  ww = Inverse[gg] - zz;
 
   uu = SparseArray[
     Flatten @ {
@@ -75,19 +67,23 @@ WickTimeReversalMoment[alpha_, {grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}] :=
 
   dgn = CirclePlus[ww[[;;dd, ;;dd]], ww[[dd+1;;, dd+1;;]]];
   off = ArrayFlatten @ {
-    {zr, ww[[;;dd, dd+1;;]]},
-    {ww[[dd+1;;, ;;dd]], zr}
+    {oo, ww[[;;dd, dd+1;;]]},
+    {ww[[dd+1;;, ;;dd]], oo}
   };
-  pf = Pfaffian[(off - zz).xx] / pf;
+  pf = Check[
+    Pfaffian[(off - zz).xx] / pf, 
+    Message[WickTimeReversalMoment::infy];
+    PseudoDivide[Pfaffian[(off - zz).xx], pf],
+    Power::infy
+  ];
 
   (* effective \Omega of \Xi *)
-  ww = Check[Inverse[zz - off], PseudoInverse[zz - off]];
+  ww = Inverse[zz - off];
   ww = off + dgn . ww . dgn;
   pf = pf * Pfaffian[xx.(ww + zz)];
 
   (* effective \Gamma *)
-  (* gg = Inverse[ww + zz]; *)
-  gg = Check[Inverse[ww + zz], PseudoInverse[ww + zz]];
+  gg = Inverse[ww + zz];
   (* effective Green's function Gij *)
   gg = CircleTimes[ThePauli[10], id] - gg;
 
@@ -96,9 +92,39 @@ WickTimeReversalMoment[alpha_, {grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}] :=
   Total[Log[2, Power[gg, alpha] + Power[1-gg, alpha]]] + Log[2, Power[pf, alpha]]
 ]
 
+
+(* NOtE: mat is supposed to be Hermitian. *)
+tameMatrix[mat_?MatrixQ, eps_] := Module[
+  {val, vec},
+  {val, vec} = Eigensystem[mat];
+  val = val /. {x_?(ZeroQ[#, eps] &) :> x + 1.25*eps};
+  Transpose[vec] . DiagonalMatrix[val] . Conjugate[vec]
+]
+
 (**** </WickTimeReversalMoment> ****)
 
 
+(**** <WickLogarithmicNegtivity> ****)
+
+(* See also: Shapourian and Ryu (2017, 2019) *)
+
+WickLogarithmicNegativity::usage = "WickLogarithmicNegativity[{grn,anm}, {k1,k2,\[Ellipsis]}] returns the fermionic negativity over modes k1, k2, \[Ellipsis] of the fermionic Gaussian states characterized by the Green's function grn (in an L\[Times]L matrix for L fermion modes) and anomalous Green's function anm (also in an L\[Times]L matrix).\nWickLogarithmicNegativity[grn, {k1,k2,\[Ellipsis]}] is equivalent to WickLogarithmicNegativity[{grn,0}, {k1,k2,\[Ellipsis]}]."
+
+Options[WickLogarithmicNegativity] = Options[WickTimeReversalMoment]
+
+
+WickLogarithmicNegativity[kk:{__Integer}][any_] :=
+  WickLogarithmicNegativity[any, kk, "Epsilon" -> OptionValue["Epsilon"]]
+
+WickLogarithmicNegativity[grn_?MatrixQ, kk:{__Integer}, opts:OptionsPattern[]] :=
+  WickTimeReversalMoment[1/2, grn, kk, opts, "Epsilon" -> OptionValue["Epsilon"]]
+
+WickLogarithmicNegativity[{grn_?MatrixQ, anm_?MatrixQ}, kk:{__Integer}, opts:OptionsPattern[]] :=
+  WickTimeReversalMoment[1/2, {grn, anm}, kk, opts, "Epsilon" -> OptionValue["Epsilon"]]
+
+(**** </WickLogarithmicNegtivity> ****)
+
+
 (**** <NambuTransforms> ****)
 
 NambuTransforms /:
@@ -360,18 +386,21 @@ WickExpectation[ws_WickState][ops:{__?AnyFermionQ}] := Module[
 
 WickGreensFunction::usage = "WickGreensFunction[ws] returns the pair {G, F} of normal Green's function G and anomalous Green's function F with respect to the Wick state ws."
 
-WickGreensFunction[ws_WickState] := Module[
+WickGreensFunction[ws_WickState] :=
+  WickGreensFunction[ws, Last @ ws]
+
+WickGreensFunction[ws_WickState, dd:{___?FermionQ}] := Module[
   { tt = ws["Stages"],
     cc = ws["Bare modes"],
     aa, bb, pp, uv, mm, ff, gg, L, n, i, j },
   {pp, uv, mm} = List @@ Most[ws];
 
-  n = Length[cc];
+  n = Length[dd];
   L = Length[pp] / 2;
 
   (* encoded modes *)
-  aa = theWickCode[cc][cc, {tt}];
-  bb = dagWickCode /@ aa;
+  aa = theWickCode[cc][dd, {tt}];
+  bb = dagWickCode /@ aa; (* NOT dagWickCode[aa] *)
 
   (* Green's function *)
   gg = Flatten @ Table[
@@ -400,17 +429,38 @@ WickGreensFunction[ws_WickState] := Module[
 WickLogarithmicNegativity[dd:{__?FermionQ}][ws_WickState] :=
   WickLogarithmicNegativity[ws, dd]
 
-WickLogarithmicNegativity[ws_WickState, dd:{__?FermionQ}] := Module[
+WickLogarithmicNegativity[ws_WickState, dd:{__?FermionQ}, opts:OptionsPattern[]] := Module[
   { gg = WickGreensFunction[ws],
     cc = Last[ws],
     kk },
   kk = dd /. Thread[cc -> Range[Length @ cc]];
-  WickTimeReversalMoment[1/2, gg, kk]
+  WickTimeReversalMoment[1/2, gg, kk, opts, "Epsilon" -> OptionValue["Epsilon"]]
 ]
 
 (**** </WickLogarithmicNegtivity> ****)
 
 
+(**** <EntanglementEntropy> ****)
+
+(* See, e.g., Peschel (2003). *)
+EntanglementEntropy[ws_WickState, dd:{__?FermionQ}] := Module[
+  { gg = WickGreensFunction[ws, dd],
+    ee = 1.25*^-5,
+    id, ff },
+  id = One[Dimensions @ First @ gg];
+  gg = ArrayFlatten @ { gg,
+    {Topple @ Last @ gg, id - Transpose[First @ gg]}
+  };
+  id = One[Dimensions @ gg];
+  ff = id - gg;
+  ff = ff . Check[ MatrixLog[ff], MatrixLog[tameMatrix[ff, ee]] ];
+  gg = gg . Check[ MatrixLog[gg], MatrixLog[tameMatrix[gg, ee]] ];
+  - (Tr[ff] + Tr[gg]) / 2
+]
+
+(**** </EntanglementEntropy> ****)
+
+
 (**** <getWickMatrix> ****)
 
 getWickMatrix::usage = "getWickMatrix[{c1,c2,...}][spec] parse spec and create the Wick matrix."
@@ -644,7 +694,7 @@ Module[
 
   (* initial state with half-filling *)
   n = Length[cc];
-  in = PadRight[Riffle[Table[1, Floor[n/2]], Table[0, Floor[n/2]]], n];
+  in = PadRight[Riffle[Table[1, Floor[n/2]], Table[0, Floor[n/2]]], n, 1];
   in = wu[WickState @ Ket[cc -> in]];
 
   (* quantum circuit with randomly selected measurements *)

Q3/PacletInfo.wl

@@ -1,6 +1,6 @@
 Paclet[
   "Name" -> "Q3",
-  "Version" -> "3.4.0",
+  "Version" -> "3.4.1",
   "WolframVersion" -> "12.3+",
   "Updating" -> Automatic,
   "Loading" -> "Startup",