Public release

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

@@ -1,6 +1,6 @@
 <|"IndexedPaths" -> None, "Driver" -> "Lucene", 
- "CreationDate" -> DateObject[{2024, 8, 23, 19, 59, 
-    50.725451`8.457800892557223}, "Instant", "Gregorian", 9.], 
+ "CreationDate" -> DateObject[{2024, 8, 24, 12, 32, 
+    31.651682`8.252971774486385}, "Instant", "Gregorian", 9.], 
  "Version" -> 5, "Synonyms" -> None, "Method" -> "BM25", 
  "Language" -> "English", "ContentFieldOptions" -> 
   <|"Title" -> <|"Stored" -> True, "Weight" -> 2|>, 

Q3/Kernel/NoisyWick.wl

@@ -1,17 +1,41 @@
 BeginPackage["Q3`"]
 
 { NoisyWickSimulate };
+{ WickDampingOperator };
 
 Begin["`Private`"]
 
+(**** <WickDampingOperator> ****)
+
+WickDampingOperator::usage = "WickDampingOperator[jmp] returns a pair {mat, const} of the quadratic kernel mat and remaining constant term const of the effective damping operator in the normal ordering that corresponds to the list jmp of quantum jump operators."
+
+WickDampingOperator::jmp = "Invalid form of quantum jump operators ``."
+
+WickDampingOperator[jmp:{__WickOperator}] :=
+  WickDampingOperator @ Flatten[ First /@ jmp ]
+
+WickDampingOperator[jmp:{Rule[_, _?VectorQ]..}] := Module[
+  { dmp = GroupBy[jmp, First, Values],
+    aa, bb },
+  If[ !ContainsOnly[Keys @ dmp, {Identity, Dagger}] || !MatrixQ[Catenate @ dmp],
+    Message[WickDampingOperator::jmp, Short @ jmp];
+    Return[$Failed]
+  ];
+
+  aa = Topple[dmp[Identity]] . dmp[Identity];
+  bb = Topple[dmp[Dagger]] . dmp[Dagger];
+  { aa - Transpose[bb], Tr[bb] } / 2
+]
+
+(**** </WickDampingOperator> ****)
+
+
 (**** <NoisyWickSimulate> ****)
 
 NoisyWickSimulate::usage = "NoisyWickSimulate[ham, jmp, in, {n, dt}] solves the quantum master equation for a non-interacting dissipative fermionic many-body system by using the Monte Carlo simulation method (alos known as the quantum jump approach or quantum trajectory method). The model is specified by the single-particle Hamiltonian matrix ham and the list jmp of quantum jump operators. The simulation starts from the initial WickState in at time 0 and goes n time steps by interval dt."
 
 NoisyWickSimulate::ham = "The Hamiltonian matrix `` needs to be numeric."
 
-NoisyWickSimulate::jmp = "Invalid form of quantum jump operators ``."
-
 NoisyWickSimulate::null = "The null state is encountered."
 
 NoisyWickSimulate::save = "The result could not be saved."
@@ -26,28 +50,21 @@ Options[NoisyWickSimulate] = {
 
 NoisyWickSimulate[ham_, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}, OptionsPattern[]] :=
   Module[
-    { dmp = GroupBy[Flatten[First /@ jmp], First, Values],
-      n = OptionValue["Samples"],
+    { n = OptionValue["Samples"],
       k = 0,
       progress = 0,
-      aa, bb, non },
-
-    If[ !ContainsOnly[Key @ dmp, {Identity, Dagger}] || !MatrixQ[Catenate @ dmp],
-      Message[NoisyWickSimulate::jmp, Short @ jmp];
-      Return[$Failed]
-    ];
+      aa, bb, dmp, fac, non },
 
-    aa = Topple[dmp[Identity]].dmp[Identity];
-    bb = Topple[dmp[Dagger]].dmp[Dagger];
-    dmp = (aa - Transpose[bb] + Tr[bb] One[Dimensions @ bb]) / 2;
+    {dmp, fac} = WickDampingOperator[jmp];
+    fac = Exp[-dt*fac];
     non = ham - I * dmp;
     non = WickGaussian @ {MatrixExp[-I*dt*non], MatrixExp[+I*dt*non]};
 
     PrintTemporary @ ProgressIndicator @ Dynamic[progress];
     data = Table[
       progress = ++k / n;
       (* theNoisyWickSimulate[non, jmp, in, {nT, dt}], *)
-      altNoisyWickSimulate[non, jmp, in, {nT, dt}],
+      altNoisyWickSimulate[{non, fac}, jmp, in, {nT, dt}],
       n
     ];
 
@@ -75,7 +92,7 @@ NoisyWickSimulate[ham_, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}, O
     False
   ]
 
-altNoisyWickSimulate[non_WickGaussian, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}] :=
+altNoisyWickSimulate[{non_WickGaussian, fac_}, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}] :=
   Module[
     { res = {in},
       new = N[in],
@@ -87,7 +104,7 @@ altNoisyWickSimulate[non_WickGaussian, jmp:{__WickOperator}, in_WickState, {nT_I
 
       (* non-unitary evolution *)
       out = non[new];
-      If[ pp < NormSquare[out],
+      If[ pp < NormSquare[fac * out],
         new = Normalize @ out;
         AppendTo[res, new];
         t += 1;
@@ -120,7 +137,7 @@ altNoisyWickSimulate[non_WickGaussian, jmp:{__WickOperator}, in_WickState, {nT_I
     Return[res]
   ]
 
-theNoisyWickSimulate[non_WickGaussian, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}] :=
+theNoisyWickSimulate[{non_WickGaussian, fac_}, jmp:{__WickOperator}, in_WickState, {nT_Integer, dt_}] :=
   Module[
     { res = {in},
       new = N[in],
@@ -132,7 +149,7 @@ theNoisyWickSimulate[non_WickGaussian, jmp:{__WickOperator}, in_WickState, {nT_I
 
       (* non-unitary evolution *)
       While[ t <= nT,
-        out = non[new];
+        out = fac * non[new];
         If[ pp < NormSquare[out],
           (* then *)
           AppendTo[res, Normalize @ out];

Q3/Kernel/Wick.wl

@@ -173,29 +173,47 @@ MakeBoxes[ws:WickState[ops:{___Rule}, cc:{___?FermionQ}], fmt_] :=
     "Interpretable" -> Automatic
   ]
 
+WickState /: (* vacuum state times a constant *)
+MakeBoxes[ws:WickState[z:Except[_?ListQ|_?ArrayQ], cc:{___?FermionQ}], fmt_] :=
+  BoxForm`ArrangeSummaryBox[
+    WickState, ws, None,
+    { BoxForm`SummaryItem @ { "Modes: ", cc },
+      BoxForm`SummaryItem @ { "Type: ", Switch[z, 0, Null, _, Vacuum] }
+    },
+    { BoxForm`SummaryItem @ { "Normalization: ", z }
+    },
+    fmt,
+    "Interpretable" -> Automatic
+  ]
+
 (* canonicalization *)
 
-WickState[cc:{__?FermionQ}] := WickState[{}, cc] (* vacuum state *)
+WickState[cc:{__?FermionQ}] := WickState[1, cc] (* vacuum state *)
+(* WickState[0, cc] represents the null state;  *)
+(* WickState[z, cc] represents a vacuum state times the overall constant z. *)
 
-WickState[0, cc:{__?FermionQ}] := (* null state *)
-  WickState[{Identity -> Table[0, Length @ cc]}, cc]
+WickState[{}, cc:{__?FermionQ}] := WickState[1, cc] (* vacuum state *)
+(* NOTE: The left-hand form cannot handle the multiplication by a global factor. *)
 
 WickState[Ket[aa_Association]] := Module[
   { cc = Select[Keys @ aa, FermionQ],
     dd = Select[Keys @ theKetTrim @ aa, FermionQ] },
   WickState[theWickOperator[Dagger @ dd, cc], cc]
 ]
 
+WickState @ Rule[cc:{__?FermionQ}, vv_/;VectorQ[vv, BinaryQ]] := Module[
+  { ww = PadRight[vv, Length @ cc, vv],
+    kk },
+  kk = PositionIndex[ww];
+  If[ MissingQ @ kk[1],
+    WickState[1, cc],
+    WickState[theWickOperator[Dagger @ cc[[kk @ 1]], cc], cc]
+  ]
+]
+
 WickState[ops:{__Dagger?AnyFermionQ}, cc:{__?FermionQ}] :=
   WickState @ Ket[cc -> 0, Peel[ops] -> 1]
 
-WickState[ops:{Rule[_, _?VectorQ]..}, cc:{___?FermionQ}] :=
-  WickState[Thread[Map[Hood, Keys @ ops] -> Values[ops]], cc] /;
-  AnyTrue[Keys @ ops, AnyFermionQ]
-
-
-WickState /:
-NormSquare[WickState[{}, cc:{___?FermionQ}]] = 1
 
 WickState /:
 NormSquare[WickState[ops:{__Rule}, cc:{__?FermionQ}]] := Quiet[
@@ -204,24 +222,45 @@ NormSquare[WickState[ops:{__Rule}, cc:{__?FermionQ}]] := Quiet[
 ]
 
 WickState /:
-Norm[ws:WickState[{___Rule}, {__?FermionQ}]] := Sqrt[NormSquare @ ws]
+NormSquare[WickState[z:Except[_?ListQ|_?ArrayQ], cc:{__?FermionQ}]] :=
+  AbsSquare[z]
 
 
 WickState /:
-Normalize[ws:WickState[{}, {__?FermionQ}], ___] = ws
+Norm[ws:WickState[{___Rule}, {__?FermionQ}]] := Sqrt[NormSquare @ ws]
+
+WickState /:
+Norm[ws:WickState[z:Except[_?ListQ|_?ArrayQ], {__?FermionQ}]] := Abs[z]
+
 
 WickState /:
 Normalize[ws:WickState[ops:{___Rule}, cc:{___?FermionQ}]] := Module[
-  { flg = Keys[ops],
+  { tag = Keys[ops],
     new },
   Quiet @ Check[
     new = Values[ops] * Power[Norm @ ws, -1/Length[ops]],
-    flg = {Identity};
+    tag = {Identity};
     new = Zero @ {1, Length @ cc}
   ];
-  WickState[Thread[flg -> new], cc]
+  WickState[Thread[tag -> new], cc]
 ]
 
+WickState /:
+Normalize[ws:WickState[z:Except[_?ListQ|_?ArrayQ], cc:{___?FermionQ}], ___] =
+  WickState[1, cc]
+
+
+WickState /:
+Times[z_, WickState[ops:{___Rule}, cc:{___?FermionQ}]] := Module[
+  { new },
+  new = Values[ops] * Power[z, 1/Length[ops]];
+  WickState[Thread[Keys[ops] -> new], cc]
+]
+
+WickState /:
+Times[z_, WickState[x:Except[_?ListQ|_?ArrayQ], cc:{___?FermionQ}]] :=
+  WickState[z x, cc]
+
 
 WickState /:
 Expand[ws:WickState[ops:{___Rule}, cc:{___?FermionQ}]] :=
@@ -237,6 +276,9 @@ Elaborate[ws:WickState[{___Rule}, cc:{__?FermionQ}]] := Module[
   KetChop[ Multiply @@ Append[ff, Ket[cc]] ]
 ]
 
+WickState /:
+Elaborate[WickState[z:Except[_?ListQ|_?ArrayQ], cc:{__?FermionQ}]] :=
+  z * Ket[cc]
 
 WickState /:
 Matrix[ws_WickState, rest___] := (
@@ -298,6 +340,11 @@ MakeBoxes[WickGaussian[{mat_?MatrixQ, inv_?MatrixQ}, rest___], fmt_] := Module[
 WickGaussian[mat_?MatrixQ, rest___] := WickGaussian[{mat, Inverse @ mat}, rest] /; 
   If[MatrixQ[mat, NumericQ], True, Message[WickGaussian::num, mat]; False]
 
+
+WickGaussian[{mat_/;MatrixQ[mat, NumericQ], inv_/;MatrixQ[inv, NumericQ]}, ___][
+  ws:WickState[Except[_?ListQ|_?ArrayQ], _]   
+] = ws
+
 WickGaussian[{mat_/;MatrixQ[mat, NumericQ], inv_/;MatrixQ[inv, NumericQ]}, ___][ws_WickState] := 
   Module[
     { new },
@@ -308,6 +355,7 @@ WickGaussian[{mat_/;MatrixQ[mat, NumericQ], inv_/;MatrixQ[inv, NumericQ]}, ___][
     WickState[new, Last @ ws]
   ]
 
+
 WickGaussian /:
 MatrixForm[op : WickGaussian[{mat_?MatrixQ, inv_?MatrixQ}, ___]] :=
   MatrixForm @ mat
@@ -455,6 +503,10 @@ Multiply[pre___, wu_WickUnitary, fs_Ket] := Multiply[pre, wu @ WickState @ fs]
 
 WickUnitary[{}, ___][any_] = any
 
+WickUnitary[mat_/;MatrixQ[mat, NumericQ], ___][
+  ws:WickState[Except[_?ListQ|_?ArrayQ], _]   
+] = ws
+
 WickUnitary[mat_/;MatrixQ[mat, NumericQ], ___][ws_WickState] := Module[
   { new },
   new = MapApply[
@@ -468,7 +520,7 @@ WickUnitary[spec__][fs_Ket] := WickUnitary[spec][WickState @ fs]
 
 
 WickUnitary /:
-ParseGate[WickUnitary[trs_, ff:{___?FermionQ}, opts___?OptionQ], more___?OptionQ] :=
+ParseGate[WickUnitary[trs_, ff:{__?FermionQ}, opts___?OptionQ], more___?OptionQ] :=
   Gate[ff, more, opts, "Label" -> "U"]
 
 (**** </WickUnitary> ****)
@@ -512,6 +564,8 @@ WickElements[ops:{___Rule}, cc:{__?FermionQ}] :=
 WickElements[WickState[ops:{___Rule}, cc:{__?FermionQ}], ___] :=
   WickElements[ops, cc]
 
+WickElements[WickState[z:Except[_?ListQ|_?ArrayQ], cc:{__?FermionQ}], ___] = {}
+
 WickElements[WickOperator[ops:{___Rule}, ___?OptionQ], cc:{__?FermionQ}] :=
   WickElements[ops, cc]
 
@@ -582,13 +636,14 @@ WickOperator[ops:{__Rule}] :=
 
 WickOperator[{}][any_] = any
 
-WickOperator[ops:{Rule[_, _?VectorQ]..}, ___][in_WickState] :=
-  WickState[Join[ops, First @ in], Last @ in]
+WickOperator[ops:{Rule[_, _?VectorQ]..}, ___][WickState[trs:{___Rule}, cc:{___?FermionQ}]] :=
+  WickState[Join[ops, trs], cc]
 
-WickOperator[ops:{__?AnyFermionQ}, ___][in_WickState] := WickState[
-  Join[theWickOperator[ops, Last @ in], First @ in],
-  Last @ in
-]
+WickOperator[ops:{Rule[_, _?VectorQ]..}, ___][WickState[z:Except[_?ListQ|_?ArrayQ], cc:{___?FermionQ}]] :=
+  z * WickState[ops, cc]
+
+WickOperator[ops:{__?AnyFermionQ}, ___][in_WickState] :=
+  WickOperator[ops, Last @ in][in]
 
 
 WickOperator /:
@@ -602,6 +657,14 @@ NonCommutativeQ[_WickOperator] = True
 WickOperator /:
 MultiplyKind[_WickOperator] = Fermion
 
+WickOperator /:
+Multiply[pre___, ops:Repeated[WickOperator[_, cc:{__?FermionQ}], {2, Infinity}], post___] :=
+  Multiply[pre, WickOperator[Flatten[First /@ {ops}], cc], post]
+
+WickOperator /:
+Multiply[pre___, aa_WickOperator, bb__WickOperator, post___] :=
+  Multiply[pre, WickOperator[Flatten[First /@ {aa, bb}]], post]
+
 WickOperator /:
 Multiply[pre___, opr_WickOperator, ws_WickState] := Multiply[pre, opr[ws]]
 
@@ -719,6 +782,11 @@ theMeasurement[ws:WickState[trs:{___Rule}, cc_], c_?FermionQ] := Module[
   ]
 ]
 
+theMeasurement[ws:WickState[z:Except[_?ListQ|_?ArrayQ], cc_], c_?FermionQ] := (
+  $MeasurementOut[c] = 0;
+  WickState[1, cc]
+)
+
 (**** </Measurement> ****)
 
 
@@ -741,13 +809,16 @@ WickExpectation[ws_][expr_Plus] :=
 WickExpectation[ws_WickState][HoldPattern @ Multiply[ops__?AnyFermionQ]] :=
   WickExpectation[ws] @ WickOperatorFrom[{ops}, Last @ ws]
 
-WickExpectation[ws:WickState[bb_, cc_]][WickOperator[ops:{___Rule}, ___]] := Module[
+WickExpectation[ws:WickState[bb:{___Rule}, cc_]][WickOperator[ops:{___Rule}, ___]] := Module[
   { aa = theConjugateReverse[bb],
     mat },
   mat = WickMatrix @ Join[aa, ops, bb];
   Pfaffian[mat] (* NOTE: The Wick state is assumed to be normalized. *)
 ]
 
+WickExpectation[ws:WickState[z:Except[_?ListQ|_?VectorQ], cc_]][op:WickOperator[ops:{___Rule}, ___]] :=
+  AbsSquare[z] * VacuumExpectation[op]
+
 (**** </WickExpectation> ****)
 
 
@@ -760,7 +831,7 @@ WickGreenFunction::null = "The null state is encountered: ``."
 WickGreenFunction[ws_WickState] :=
   WickGreenFunction[ws, Last @ ws]
 
-WickGreenFunction[ws:WickState[qq_, cc_], dd:{___?FermionQ}] := Module[
+WickGreenFunction[ws:WickState[qq:{__Rule}, cc_], dd:{___?FermionQ}] := Module[
   { pp = theConjugateReverse[qq],
     aa, bb, gg, wm, n },
   bb = Lookup[First /@ PositionIndex[cc], dd];
@@ -801,6 +872,9 @@ WickGreenFunction[ws:WickState[qq_, cc_], dd:{___?FermionQ}] := Module[
   Return[gg] (* NOTE: The Wick state is assumed to be normalized. *)
 ]
 
+WickGreenFunction[WickState[z:Except[_?ListQ|_?ArrayQ], cc_], dd:{___?FermionQ}] :=
+  One[Length @ dd]
+
 (**** </WickGreenFunction> ****)
 
 

Q3/PacletInfo.wl

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

RELEASES.md

@@ -1,5 +1,10 @@
 # Release Notes
 
+## 3.5.8
+
+- New functions: WickDampingOperator
+- Bug fixes: NoisyWickSimulate
+
 ## 3.5.7
 
 - Improved: WickGreenFunction