diff --git a/tensorflow_quantum/core/src/BUILD b/tensorflow_quantum/core/src/BUILD
index 94dffccf9..65a5edb4f 100644
--- a/tensorflow_quantum/core/src/BUILD
+++ b/tensorflow_quantum/core/src/BUILD
@@ -15,6 +15,34 @@ cc_library(
     ],
 )
 
+cc_library(
+    name = "adj_hessian_util",
+    srcs = ["adj_hessian_util.cc"],
+    hdrs = ["adj_hessian_util.h", "adj_util.h"],
+    deps = [
+        ":circuit_parser_qsim",
+        "@qsim//lib:circuit",
+        "@qsim//lib:fuser",
+        "@qsim//lib:fuser_basic",
+        "@qsim//lib:gate",
+        "@qsim//lib:gates_cirq",
+        "@qsim//lib:io",
+        "@qsim//lib:matrix",
+    ],
+)
+
+cc_test(
+    name = "adj_hessian_util_test",
+    srcs = ["adj_hessian_util_test.cc"],
+    deps = [
+        ":adj_hessian_util",
+        ":circuit_parser_qsim",
+        "@com_google_googletest//:gtest_main",
+        "@qsim//lib:gates_cirq",
+        "@qsim//lib:matrix",
+    ],
+)
+
 cc_library(
     name = "adj_util",
     srcs = ["adj_util.cc"],
diff --git a/tensorflow_quantum/core/src/adj_hessian_util.cc b/tensorflow_quantum/core/src/adj_hessian_util.cc
new file mode 100644
index 000000000..8b4b122a4
--- /dev/null
+++ b/tensorflow_quantum/core/src/adj_hessian_util.cc
@@ -0,0 +1,482 @@
+/* Copyright 2021 The TensorFlow Quantum Authors. All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+==============================================================================*/
+#include "tensorflow_quantum/core/src/adj_hessian_util.h"
+
+#include <functional>
+#include <string>
+#include <vector>
+
+#include "../qsim/lib/circuit.h"
+#include "../qsim/lib/fuser.h"
+#include "../qsim/lib/fuser_basic.h"
+#include "../qsim/lib/gate.h"
+#include "../qsim/lib/gates_cirq.h"
+#include "../qsim/lib/io.h"
+#include "../qsim/lib/matrix.h"
+#include "tensorflow_quantum/core/src/circuit_parser_qsim.h"
+
+namespace tfq {
+
+typedef qsim::Cirq::GateCirq<float> QsimGate;
+typedef qsim::Circuit<QsimGate> QsimCircuit;
+
+void CreateHessianCircuit(
+    const QsimCircuit& circuit, const std::vector<GateMetaData>& metadata,
+    std::vector<std::vector<qsim::GateFused<QsimGate>>>* partial_fuses,
+    std::vector<GradientOfGate>* grad_gates) {
+  for (std::vector<tfq::GateMetaData>::size_type i = 0; i < metadata.size();
+       i++) {
+    if (metadata[i].symbol_values.empty()) {
+      continue;
+    }
+    // found a gate that was constructed with symbols.
+    GradientOfGate grad;
+
+    // Single qubit Eigen.
+    if (circuit.gates[i].kind == qsim::Cirq::GateKind::kXPowGate ||
+        circuit.gates[i].kind == qsim::Cirq::GateKind::kYPowGate ||
+        circuit.gates[i].kind == qsim::Cirq::GateKind::kZPowGate ||
+        circuit.gates[i].kind == qsim::Cirq::GateKind::kHPowGate) {
+      PopulateHessianSingleEigen(
+          metadata[i].create_f1, metadata[i].symbol_values[0], i,
+          circuit.gates[i].qubits[0], metadata[i].gate_params[0],
+          metadata[i].gate_params[1], metadata[i].gate_params[2], &grad);
+      grad_gates->push_back(grad);
+    }
+
+    // Two qubit Eigen.
+    else if (circuit.gates[i].kind == qsim::Cirq::GateKind::kCZPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kCXPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kXXPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kYYPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kZZPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kISwapPowGate ||
+             circuit.gates[i].kind == qsim::Cirq::GateKind::kSwapPowGate) {
+      bool swapq = circuit.gates[i].swapped;
+      PopulateHessianTwoEigen(
+          metadata[i].create_f2, metadata[i].symbol_values[0], i,
+          swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+          swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+          metadata[i].gate_params[0], metadata[i].gate_params[1],
+          metadata[i].gate_params[2], &grad);
+      grad_gates->push_back(grad);
+    }
+
+    // PhasedX
+    else if (circuit.gates[i].kind == qsim::Cirq::GateKind::kPhasedXPowGate) {
+      // Process potentially several symbols.
+      bool symbolic_pexp = false;
+      bool symbolic_exp = false;
+      for (std::vector<std::basic_string<char>>::size_type j = 0;
+           j < metadata[i].symbol_values.size(); j++) {
+        if (metadata[i].placeholder_names[j] ==
+            GateParamNames::kPhaseExponent) {
+          symbolic_pexp = true;
+          PopulateHessianPhasedXPhasedExponent(
+              metadata[i].symbol_values[j], i, circuit.gates[i].qubits[0],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3],
+              metadata[i].gate_params[4], &grad);
+        } else if (metadata[i].placeholder_names[j] ==
+                   GateParamNames::kExponent) {
+          symbolic_exp = true;
+          PopulateHessianPhasedXExponent(
+              metadata[i].symbol_values[j], i, circuit.gates[i].qubits[0],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3],
+              metadata[i].gate_params[4], &grad);
+        }
+      }
+      if (symbolic_pexp && symbolic_exp) {
+        PopulateCrossTermPhasedXPhasedExponentExponent(
+            i, circuit.gates[i].qubits[0], metadata[i].gate_params[0],
+            metadata[i].gate_params[1], metadata[i].gate_params[2],
+            metadata[i].gate_params[3], metadata[i].gate_params[4], &grad);
+      }
+      grad_gates->push_back(grad);
+    }
+
+    // Fsim
+    else if (circuit.gates[i].kind == qsim::Cirq::GateKind::kFSimGate) {
+      // Process potentially several symbols.
+
+      bool swapq = circuit.gates[i].swapped;
+      bool symbolic_theta = false;
+      bool symbolic_phi = false;
+      for (std::vector<std::basic_string<char>>::size_type j = 0;
+           j < metadata[i].symbol_values.size(); j++) {
+        if (metadata[i].placeholder_names[j] == GateParamNames::kTheta) {
+          symbolic_theta = true;
+          PopulateHessianFsimTheta(
+              metadata[i].symbol_values[j], i,
+              swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+              swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+        } else if (metadata[i].placeholder_names[j] == GateParamNames::kPhi) {
+          symbolic_phi = true;
+          PopulateHessianFsimPhi(
+              metadata[i].symbol_values[j], i,
+              swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+              swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+        }
+      }
+      if (symbolic_theta && symbolic_phi) {
+        PopulateCrossTermFsimThetaPhi(
+            i, swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+            swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+            metadata[i].gate_params[0], metadata[i].gate_params[1],
+            metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+      }
+
+      grad_gates->push_back(grad);
+    }
+
+    // PhasedISwap
+    else if (circuit.gates[i].kind ==
+             qsim::Cirq::GateKind::kPhasedISwapPowGate) {
+      // Process potentially several symbols.
+      bool swapq = circuit.gates[i].swapped;
+      bool symbolic_pexp = false;
+      bool symbolic_exp = false;
+      for (std::vector<std::basic_string<char>>::size_type j = 0;
+           j < metadata[i].symbol_values.size(); j++) {
+        if (metadata[i].placeholder_names[j] ==
+            GateParamNames::kPhaseExponent) {
+          symbolic_pexp = true;
+          PopulateHessianPhasedISwapPhasedExponent(
+              metadata[i].symbol_values[j], i,
+              swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+              swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+
+        } else if (metadata[i].placeholder_names[j] ==
+                   GateParamNames::kExponent) {
+          symbolic_exp = true;
+          PopulateHessianPhasedISwapExponent(
+              metadata[i].symbol_values[j], i,
+              swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+              swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+              metadata[i].gate_params[0], metadata[i].gate_params[1],
+              metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+        }
+      }
+      if (symbolic_pexp && symbolic_exp) {
+        PopulateCrossTermPhasedISwapPhasedExponentExponent(
+            i, swapq ? circuit.gates[i].qubits[1] : circuit.gates[i].qubits[0],
+            swapq ? circuit.gates[i].qubits[0] : circuit.gates[i].qubits[1],
+            metadata[i].gate_params[0], metadata[i].gate_params[1],
+            metadata[i].gate_params[2], metadata[i].gate_params[3], &grad);
+      }
+      grad_gates->push_back(grad);
+    }
+  }
+
+  // Produce partial fuses around the hessian gates.
+  auto fuser = qsim::BasicGateFuser<qsim::IO, QsimGate>();
+  auto left = circuit.gates.begin();
+  auto right = left;
+
+  partial_fuses->assign(grad_gates->size() + 1,
+                        std::vector<qsim::GateFused<QsimGate>>({}));
+  for (std::vector<GradientOfGate>::size_type i = 0; i < grad_gates->size();
+       i++) {
+    right = circuit.gates.begin() + (*grad_gates)[i].index;
+    (*partial_fuses)[i] =
+        fuser.FuseGates(qsim::BasicGateFuser<qsim::IO, QsimGate>::Parameter(),
+                        circuit.num_qubits, left, right);
+    left = right + 1;
+  }
+  right = circuit.gates.end();
+  (*partial_fuses)[grad_gates->size()] =
+      fuser.FuseGates(qsim::BasicGateFuser<qsim::IO, QsimGate>::Parameter(),
+                      circuit.num_qubits, left, right);
+}
+
+void PopulateHessianSingleEigen(
+    const std::function<QsimGate(unsigned int, unsigned int, float, float)>&
+        create_f,
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    float exp, float exp_s, float gs, GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = create_f(0, qid, (exp + _HESS_EPS) * exp_s, gs);
+  auto center = create_f(0, qid, exp * exp_s, gs);
+  auto right = create_f(0, qid, (exp - _HESS_EPS) * exp_s, gs);
+  // Due to precision issue, (1) multiplies weights first rather than last.
+  // and (2) doesn't use _INVERSE_HESS_EPS_SQUARE
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix2Add(right.matrix,
+             left.matrix);  // left's entries have right added
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix2Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianTwoEigen(
+    const std::function<QsimGate(unsigned int, unsigned int, unsigned int,
+                                 float, float)>& create_f,
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    unsigned int qid2, float exp, float exp_s, float gs, GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = create_f(0, qid, qid2, (exp + _HESS_EPS) * exp_s, gs);
+  auto center = create_f(0, qid, qid2, exp * exp_s, gs);
+  auto right = create_f(0, qid, qid2, (exp - _HESS_EPS) * exp_s, gs);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix4Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianPhasedXPhasedExponent(const std::string& symbol,
+                                          unsigned int location,
+                                          unsigned int qid, float pexp,
+                                          float pexp_s, float exp, float exp_s,
+                                          float gs, GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp + _HESS_EPS) * pexp_s, exp * exp_s, gs);
+  auto center = qsim::Cirq::PhasedXPowGate<float>::Create(0, qid, pexp * pexp_s,
+                                                          exp * exp_s, gs);
+  auto right = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp - _HESS_EPS) * pexp_s, exp * exp_s, gs);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix2Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix2Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianPhasedXExponent(const std::string& symbol,
+                                    unsigned int location, unsigned int qid,
+                                    float pexp, float pexp_s, float exp,
+                                    float exp_s, float gs,
+                                    GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, pexp * pexp_s, (exp + _HESS_EPS) * exp_s, gs);
+  auto center = qsim::Cirq::PhasedXPowGate<float>::Create(0, qid, pexp * pexp_s,
+                                                          exp * exp_s, gs);
+  auto right = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, pexp * pexp_s, (exp - _HESS_EPS) * exp_s, gs);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix2Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix2Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateCrossTermPhasedXPhasedExponentExponent(
+    unsigned int location, unsigned int qid, float pexp, float pexp_s,
+    float exp, float exp_s, float gs, GradientOfGate* grad) {
+  grad->params.push_back(kUsePrevTwoSymbols);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp + _GRAD_EPS) * pexp_s, (exp + _GRAD_EPS) * exp_s, gs);
+  auto left_center = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp + _GRAD_EPS) * pexp_s, (exp - _GRAD_EPS) * exp_s, gs);
+  auto right_center = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp - _GRAD_EPS) * pexp_s, (exp + _GRAD_EPS) * exp_s, gs);
+  auto right = qsim::Cirq::PhasedXPowGate<float>::Create(
+      0, qid, (pexp - _GRAD_EPS) * pexp_s, (exp - _GRAD_EPS) * exp_s, gs);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left_center.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right_center.matrix);
+  Matrix2Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  Matrix2Add(right_center.matrix, left_center.matrix);
+  Matrix2Diff(left_center.matrix,
+              left.matrix);  // left's entries have left_center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianFsimTheta(const std::string& symbol, unsigned int location,
+                              unsigned int qid, unsigned qid2, float theta,
+                              float theta_s, float phi, float phi_s,
+                              GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta + _HESS_EPS) * theta_s, phi * phi_s);
+  auto center = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, theta * theta_s, phi * phi_s);
+  auto right = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta - _HESS_EPS) * theta_s, phi * phi_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix4Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianFsimPhi(const std::string& symbol, unsigned int location,
+                            unsigned int qid, unsigned qid2, float theta,
+                            float theta_s, float phi, float phi_s,
+                            GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::FSimGate<float>::Create(0, qid, qid2, theta * theta_s,
+                                                  (phi + _HESS_EPS) * phi_s);
+  auto center = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, theta * theta_s, phi * phi_s);
+  auto right = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, theta * theta_s, (phi - _HESS_EPS) * phi_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix4Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateCrossTermFsimThetaPhi(unsigned int location, unsigned int qid,
+                                   unsigned qid2, float theta, float theta_s,
+                                   float phi, float phi_s,
+                                   GradientOfGate* grad) {
+  grad->params.push_back(kUsePrevTwoSymbols);
+  grad->index = location;
+  auto left = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta + _GRAD_EPS) * theta_s, (phi + _GRAD_EPS) * phi_s);
+  auto left_center = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta + _GRAD_EPS) * theta_s, (phi - _GRAD_EPS) * phi_s);
+  auto right_center = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta - _GRAD_EPS) * theta_s, (phi + _GRAD_EPS) * phi_s);
+  auto right = qsim::Cirq::FSimGate<float>::Create(
+      0, qid, qid2, (theta - _GRAD_EPS) * theta_s, (phi - _GRAD_EPS) * phi_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left_center.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right_center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  Matrix4Add(right_center.matrix, left_center.matrix);
+  Matrix4Diff(left_center.matrix,
+              left.matrix);  // left's entries have left_center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianPhasedISwapPhasedExponent(
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    unsigned int qid2, float pexp, float pexp_s, float exp, float exp_s,
+    GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp + _HESS_EPS) * pexp_s, exp * exp_s);
+  auto center = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, pexp * pexp_s, exp * exp_s);
+  auto right = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp - _HESS_EPS) * pexp_s, exp * exp_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix4Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateHessianPhasedISwapExponent(const std::string& symbol,
+                                        unsigned int location, unsigned int qid,
+                                        unsigned int qid2, float pexp,
+                                        float pexp_s, float exp, float exp_s,
+                                        GradientOfGate* grad) {
+  grad->params.push_back(symbol);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, pexp * pexp_s, (exp + _HESS_EPS) * exp_s);
+  auto center = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, pexp * pexp_s, exp * exp_s);
+  auto right = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, pexp * pexp_s, (exp - _HESS_EPS) * exp_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  qsim::MatrixScalarMultiply(2.0, center.matrix);
+  Matrix4Diff(center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+void PopulateCrossTermPhasedISwapPhasedExponentExponent(
+    unsigned int location, unsigned int qid, unsigned int qid2, float pexp,
+    float pexp_s, float exp, float exp_s, GradientOfGate* grad) {
+  grad->params.push_back(kUsePrevTwoSymbols);
+  grad->index = location;
+  auto left = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp + _GRAD_EPS) * pexp_s, (exp + _GRAD_EPS) * exp_s);
+  auto left_center = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp + _GRAD_EPS) * pexp_s, (exp - _GRAD_EPS) * exp_s);
+  auto right_center = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp - _GRAD_EPS) * pexp_s, (exp + _GRAD_EPS) * exp_s);
+  auto right = qsim::Cirq::PhasedISwapPowGate<float>::Create(
+      0, qid, qid2, (pexp - _GRAD_EPS) * pexp_s, (exp - _GRAD_EPS) * exp_s);
+  // Due to precision issue, multiply weights first.
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, left_center.matrix);
+  qsim::MatrixScalarMultiply(_INVERSE_HESS_EPS_SQUARE, right_center.matrix);
+  Matrix4Add(right.matrix,
+             left.matrix);  // left's entries have right added.
+  Matrix4Add(right_center.matrix, left_center.matrix);
+  Matrix4Diff(left_center.matrix,
+              left.matrix);  // left's entries have center subtracted.
+  grad->grad_gates.push_back(left);
+}
+
+}  // namespace tfq
diff --git a/tensorflow_quantum/core/src/adj_hessian_util.h b/tensorflow_quantum/core/src/adj_hessian_util.h
new file mode 100644
index 000000000..9178e071e
--- /dev/null
+++ b/tensorflow_quantum/core/src/adj_hessian_util.h
@@ -0,0 +1,126 @@
+/* Copyright 2021 The TensorFlow Quantum Authors. All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+==============================================================================*/
+
+#ifndef TFQ_CORE_SRC_ADJ_HESSIAN_UTIL_H_
+#define TFQ_CORE_SRC_ADJ_HESSIAN_UTIL_H_
+
+#include <functional>
+#include <string>
+#include <vector>
+
+#include "../qsim/lib/circuit.h"
+#include "../qsim/lib/fuser.h"
+#include "../qsim/lib/fuser_basic.h"
+#include "../qsim/lib/gates_cirq.h"
+#include "../qsim/lib/io.h"
+#include "../qsim/lib/matrix.h"
+#include "tensorflow_quantum/core/src/adj_util.h"
+#include "tensorflow_quantum/core/src/circuit_parser_qsim.h"
+
+namespace tfq {
+
+static const float _GRAD_EPS = 5e-3;
+static const float _HESS_EPS = 1e-2;
+static const float _INVERSE_HESS_EPS_SQUARE = 1e4;
+static const std::string kUsePrevTwoSymbols = "use_prev_two_symbols";
+
+// Computes all gates who's hessian will need to be taken, in addition
+// fuses all gates around those gates for faster circuit execution.
+void CreateHessianCircuit(
+    const qsim::Circuit<qsim::Cirq::GateCirq<float>>& circuit,
+    const std::vector<GateMetaData>& metadata,
+    std::vector<std::vector<qsim::GateFused<qsim::Cirq::GateCirq<float>>>>*
+        partial_fuses,
+    std::vector<GradientOfGate>* grad_gates);
+
+void PopulateHessianSingleEigen(
+    const std::function<qsim::Cirq::GateCirq<float>(unsigned int, unsigned int,
+                                                    float, float)>& create_f,
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    float exp, float exp_s, float gs, GradientOfGate* grad);
+
+void PopulateHessianTwoEigen(
+    const std::function<qsim::Cirq::GateCirq<float>(
+        unsigned int, unsigned int, unsigned int, float, float)>& create_f,
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    unsigned int qid2, float exp, float exp_s, float gs, GradientOfGate* grad);
+
+// Note: all methods below expect gate qubit indices to have been swapped so
+// qid < qid2.
+void PopulateHessianPhasedXPhasedExponent(const std::string& symbol,
+                                          unsigned int location,
+                                          unsigned int qid, float pexp,
+                                          float pexp_s, float exp, float exp_s,
+                                          float gs, GradientOfGate* grad);
+
+void PopulateHessianPhasedXExponent(const std::string& symbol,
+                                    unsigned int location, unsigned int qid,
+                                    float pexp, float pexp_s, float exp,
+                                    float exp_s, float gs,
+                                    GradientOfGate* grad);
+
+void PopulateCrossTermPhasedXPhasedExponentExponent(
+    unsigned int location, unsigned int qid, float pexp, float pexp_s,
+    float exp, float exp_s, float gs, GradientOfGate* grad);
+
+void PopulateHessianFsimTheta(const std::string& symbol, unsigned int location,
+                              unsigned int qid, unsigned qid2, float theta,
+                              float theta_s, float phi, float phi_s,
+                              GradientOfGate* grad);
+
+void PopulateHessianFsimPhi(const std::string& symbol, unsigned int location,
+                            unsigned int qid, unsigned qid2, float theta,
+                            float theta_s, float phi, float phi_s,
+                            GradientOfGate* grad);
+
+void PopulateCrossTermFsimThetaPhi(unsigned int location, unsigned int qid,
+                                   unsigned qid2, float theta, float theta_s,
+                                   float phi, float phi_s,
+                                   GradientOfGate* grad);
+
+void PopulateHessianPhasedISwapPhasedExponent(
+    const std::string& symbol, unsigned int location, unsigned int qid,
+    unsigned int qid2, float pexp, float pexp_s, float exp, float exp_s,
+    GradientOfGate* grad);
+
+void PopulateHessianPhasedISwapExponent(const std::string& symbol,
+                                        unsigned int location, unsigned int qid,
+                                        unsigned int qid2, float pexp,
+                                        float pexp_s, float exp, float exp_s,
+                                        GradientOfGate* grad);
+
+void PopulateCrossTermPhasedISwapPhasedExponentExponent(
+    unsigned int location, unsigned int qid, unsigned int qid2, float pexp,
+    float pexp_s, float exp, float exp_s, GradientOfGate* grad);
+
+// does matrix elementwise addition dest += source.
+template <typename Array2>
+void Matrix2Add(Array2& source, Array2& dest) {
+  for (unsigned i = 0; i < 8; i++) {
+    dest[i] += source[i];
+  }
+}
+
+// does matrix elementwise addition dest += source.
+template <typename Array2>
+void Matrix4Add(Array2& source, Array2& dest) {
+  for (unsigned i = 0; i < 32; i++) {
+    dest[i] += source[i];
+  }
+}
+
+}  // namespace tfq
+
+#endif  // TFQ_CORE_SRC_ADJ_UTIL_H_
diff --git a/tensorflow_quantum/core/src/adj_hessian_util_test.cc b/tensorflow_quantum/core/src/adj_hessian_util_test.cc
new file mode 100644
index 000000000..62c3dba87
--- /dev/null
+++ b/tensorflow_quantum/core/src/adj_hessian_util_test.cc
@@ -0,0 +1,710 @@
+/* Copyright 2020 The TensorFlow Quantum Authors. All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+==============================================================================*/
+#include "tensorflow_quantum/core/src/adj_hessian_util.h"
+
+#include <string>
+#include <vector>
+
+#include "../qsim/lib/circuit.h"
+#include "../qsim/lib/fuser.h"
+#include "../qsim/lib/fuser_basic.h"
+#include "../qsim/lib/gates_cirq.h"
+#include "../qsim/lib/io.h"
+#include "../qsim/lib/matrix.h"
+#include "gtest/gtest.h"
+
+namespace tfq {
+namespace {
+
+void Matrix2Equal(const std::vector<float>& v,
+                  const std::vector<float>& expected, float eps) {
+  for (int i = 0; i < 8; i++) {
+    EXPECT_NEAR(v[i], expected[i], eps);
+  }
+}
+
+void Matrix4Equal(const std::vector<float>& v,
+                  const std::vector<float>& expected, float eps) {
+  for (int i = 0; i < 32; i++) {
+    EXPECT_NEAR(v[i], expected[i], eps);
+  }
+}
+
+typedef absl::flat_hash_map<std::string, std::pair<int, float>> SymbolMap;
+typedef qsim::Cirq::GateCirq<float> QsimGate;
+typedef qsim::Circuit<QsimGate> QsimCircuit;
+
+class OneQubitEigenFixture
+    : public ::testing::TestWithParam<
+          std::function<QsimGate(unsigned int, unsigned int, float, float)>> {};
+
+TEST_P(OneQubitEigenFixture, CreateHessianSingleEigen) {
+  QsimCircuit circuit;
+  std::vector<GateMetaData> metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  // Create a symbolized gate.
+  std::function<QsimGate(unsigned int, unsigned int, float, float)> given_f =
+      GetParam();
+
+  circuit.num_qubits = 2;
+  circuit.gates.push_back(given_f(0, 1, 1.0, 2.0));
+  GateMetaData meta;
+  meta.index = 0;
+  meta.symbol_values.push_back("TheSymbol");
+  meta.placeholder_names.push_back(GateParamNames::kExponent);
+  meta.gate_params = {1.0, 1.0, 2.0};
+  meta.create_f1 = given_f;
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 1);
+  EXPECT_EQ(grad_gates[0].index, 0);
+  EXPECT_EQ(grad_gates[0].params.size(), 1);
+  EXPECT_EQ(grad_gates[0].params[0], "TheSymbol");
+  // fuse everything into 2 gates. One fuse before this gate and one after.
+  // both wind up being identity since this is the only gate.
+  EXPECT_EQ(fuses.size(), 2);
+
+  GradientOfGate tmp;
+  PopulateHessianSingleEigen(given_f, "TheSymbol", 0, 1, 1.0, 1.0, 2.0, &tmp);
+
+  Matrix2Equal(tmp.grad_gates[0].matrix, grad_gates[0].grad_gates[0].matrix,
+               1e-4);
+
+  // Test with NO symbol.
+  metadata.clear();
+  meta.symbol_values.clear();
+  meta.placeholder_names.clear();
+  fuses.clear();
+  grad_gates.clear();
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 0);
+  EXPECT_EQ(fuses.size(), 1);
+}
+
+INSTANTIATE_TEST_CASE_P(
+    OneQubitEigenTests, OneQubitEigenFixture,
+    ::testing::Values(&qsim::Cirq::XPowGate<float>::Create,
+                      &qsim::Cirq::YPowGate<float>::Create,
+                      &qsim::Cirq::ZPowGate<float>::Create,
+                      &qsim::Cirq::HPowGate<float>::Create));
+
+class TwoQubitEigenFixture
+    : public ::testing::TestWithParam<std::function<QsimGate(
+          unsigned int, unsigned int, unsigned int, float, float)>> {};
+
+TEST_P(TwoQubitEigenFixture, CreateHessianTwoEigen) {
+  QsimCircuit circuit;
+  std::vector<GateMetaData> metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  // Create a symbolized gate.
+  std::function<QsimGate(unsigned int, unsigned int, unsigned int, float,
+                         float)>
+      given_f = GetParam();
+
+  circuit.num_qubits = 2;
+  circuit.gates.push_back(given_f(0, 0, 1, 1.0, 2.0));
+  GateMetaData meta;
+  meta.index = 0;
+  meta.symbol_values.push_back("TheSymbol");
+  meta.placeholder_names.push_back(GateParamNames::kExponent);
+  meta.gate_params = {1.0, 1.0, 2.0};
+  meta.create_f2 = given_f;
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 1);
+  EXPECT_EQ(grad_gates[0].index, 0);
+  EXPECT_EQ(grad_gates[0].params.size(), 1);
+  EXPECT_EQ(grad_gates[0].params[0], "TheSymbol");
+
+  // fuse everything into 2 gates. One fuse before this gate and one after.
+  // both wind up being identity since this is the only gate.
+  EXPECT_EQ(fuses.size(), 2);
+
+  GradientOfGate tmp;
+  PopulateHessianTwoEigen(given_f, "TheSymbol", 0, 0, 1, 1.0, 1.0, 2.0, &tmp);
+
+  Matrix4Equal(tmp.grad_gates[0].matrix, grad_gates[0].grad_gates[0].matrix,
+               1e-4);
+
+  // Test with NO symbol.
+  metadata.clear();
+  meta.symbol_values.clear();
+  meta.placeholder_names.clear();
+  fuses.clear();
+  grad_gates.clear();
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 0);
+  EXPECT_EQ(fuses.size(), 1);  // fuse everything into 1 gate.
+}
+
+INSTANTIATE_TEST_CASE_P(
+    TwoQubitEigenTests, TwoQubitEigenFixture,
+    ::testing::Values(&qsim::Cirq::CZPowGate<float>::Create,
+                      &qsim::Cirq::CXPowGate<float>::Create,
+                      &qsim::Cirq::XXPowGate<float>::Create,
+                      &qsim::Cirq::YYPowGate<float>::Create,
+                      &qsim::Cirq::ZZPowGate<float>::Create,
+                      &qsim::Cirq::ISwapPowGate<float>::Create,
+                      &qsim::Cirq::SwapPowGate<float>::Create));
+
+TEST(AdjHessianUtilTest, CreateHessianPhasedX) {
+  QsimCircuit circuit;
+  std::vector<GateMetaData> metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  // Create a symbolized gate.
+  circuit.num_qubits = 2;
+  circuit.gates.push_back(
+      qsim::Cirq::PhasedXPowGate<float>::Create(0, 0, 1.0, 2.0, 3.0));
+  GateMetaData meta;
+  meta.index = 0;
+  meta.symbol_values.push_back("TheSymbol");
+  meta.placeholder_names.push_back(GateParamNames::kPhaseExponent);
+  meta.symbol_values.push_back("TheSymbol2");
+  meta.placeholder_names.push_back(GateParamNames::kExponent);
+  meta.gate_params = {1.0, 1.0, 2.0, 1.0, 3.0};
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 1);
+  EXPECT_EQ(grad_gates[0].index, 0);
+  EXPECT_EQ(grad_gates[0].params.size(), 3);
+  EXPECT_EQ(grad_gates[0].params[0], "TheSymbol");
+  EXPECT_EQ(grad_gates[0].params[1], "TheSymbol2");
+  // Third symbol is automatically generated `kUsePrevTwoSymbols`.
+  EXPECT_EQ(grad_gates[0].params[2], kUsePrevTwoSymbols);
+
+  // fuse everything into 2 gates. One fuse before this gate and one after.
+  // both wind up being identity since this is the only gate.
+  EXPECT_EQ(fuses.size(), 2);
+
+  GradientOfGate tmp;
+  PopulateHessianPhasedXPhasedExponent("TheSymbol", 0, 0, 1.0, 1.0, 2.0, 1.0,
+                                       3.0, &tmp);
+
+  Matrix2Equal(tmp.grad_gates[0].matrix, grad_gates[0].grad_gates[0].matrix,
+               1e-4);
+
+  GradientOfGate tmp2;
+  PopulateHessianPhasedXExponent("TheSymbol2", 0, 0, 1.0, 1.0, 2.0, 1.0, 3.0,
+                                 &tmp2);
+
+  Matrix2Equal(tmp2.grad_gates[0].matrix, grad_gates[0].grad_gates[1].matrix,
+               1e-4);
+
+  GradientOfGate tmp3;
+  PopulateHessianPhasedXExponent("TheSymbol31", 0, 0, 1.0, 1.0, 2.0, 1.0, 3.0,
+                                 &tmp2);
+
+  Matrix2Equal(tmp2.grad_gates[0].matrix, grad_gates[0].grad_gates[1].matrix,
+               1e-4);
+  metadata.clear();
+  meta.symbol_values.clear();
+  meta.placeholder_names.clear();
+  grad_gates.clear();
+  fuses.clear();
+
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 0);
+  EXPECT_EQ(fuses.size(), 1);
+}
+
+TEST(AdjHessianUtilTest, CreateHessianPhasedISwap) {
+  QsimCircuit circuit;
+  std::vector<GateMetaData> metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  // Create a symbolized gate.
+  circuit.num_qubits = 2;
+  circuit.gates.push_back(
+      qsim::Cirq::PhasedISwapPowGate<float>::Create(0, 0, 1, 1.0, 2.0));
+  GateMetaData meta;
+  meta.index = 0;
+  meta.symbol_values.push_back("TheSymbol");
+  meta.placeholder_names.push_back(GateParamNames::kPhaseExponent);
+  meta.symbol_values.push_back("TheSymbol2");
+  meta.placeholder_names.push_back(GateParamNames::kExponent);
+  meta.gate_params = {1.0, 1.0, 2.0, 1.0};
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 1);
+  EXPECT_EQ(grad_gates[0].index, 0);
+  EXPECT_EQ(grad_gates[0].params.size(), 3);
+  EXPECT_EQ(grad_gates[0].params[0], "TheSymbol");
+  EXPECT_EQ(grad_gates[0].params[1], "TheSymbol2");
+  // Third symbol is automatically generated `kUsePrevTwoSymbols`.
+  EXPECT_EQ(grad_gates[0].params[2], kUsePrevTwoSymbols);
+
+  // fuse everything into 2 gates. One fuse before this gate and one after.
+  // both wind up being identity since this is the only gate.
+  EXPECT_EQ(fuses.size(), 2);
+
+  GradientOfGate tmp;
+  PopulateHessianPhasedISwapPhasedExponent("TheSymbol", 0, 0, 1, 1.0, 1.0, 2.0,
+                                           1.0, &tmp);
+
+  Matrix4Equal(tmp.grad_gates[0].matrix, grad_gates[0].grad_gates[0].matrix,
+               3e-2);
+
+  GradientOfGate tmp2;
+  PopulateHessianPhasedISwapExponent("TheSymbol2", 0, 0, 1, 1.0, 1.0, 2.0, 1.0,
+                                     &tmp2);
+
+  Matrix4Equal(tmp2.grad_gates[0].matrix, grad_gates[0].grad_gates[1].matrix,
+               3e-2);
+
+  metadata.clear();
+  meta.symbol_values.clear();
+  meta.placeholder_names.clear();
+  grad_gates.clear();
+  fuses.clear();
+
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 0);
+  EXPECT_EQ(fuses.size(), 1);
+}
+
+TEST(AdjHessianUtilTest, CreateHessianFSim) {
+  QsimCircuit circuit;
+  std::vector<GateMetaData> metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  // Create a symbolized gate.
+  circuit.num_qubits = 2;
+  circuit.gates.push_back(
+      qsim::Cirq::FSimGate<float>::Create(0, 0, 1, 1.0, 2.0));
+  GateMetaData meta;
+  meta.index = 0;
+  meta.symbol_values.push_back("TheSymbol");
+  meta.placeholder_names.push_back(GateParamNames::kTheta);
+  meta.symbol_values.push_back("TheSymbol2");
+  meta.placeholder_names.push_back(GateParamNames::kPhi);
+  meta.gate_params = {1.0, 1.0, 2.0, 1.0};
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 1);
+  EXPECT_EQ(grad_gates[0].index, 0);
+  EXPECT_EQ(grad_gates[0].params.size(), 3);
+  EXPECT_EQ(grad_gates[0].params[0], "TheSymbol");
+  EXPECT_EQ(grad_gates[0].params[1], "TheSymbol2");
+  // Third symbol is automatically generated `kUsePrevTwoSymbols`.
+  EXPECT_EQ(grad_gates[0].params[2], kUsePrevTwoSymbols);
+
+  // fuse everything into 2 gates. One fuse before this gate and one after.
+  // both wind up being identity since this is the only gate.
+  EXPECT_EQ(fuses.size(), 2);
+
+  GradientOfGate tmp;
+  PopulateHessianFsimTheta("TheSymbol", 0, 0, 1, 1.0, 1.0, 2.0, 1.0, &tmp);
+
+  Matrix4Equal(tmp.grad_gates[0].matrix, grad_gates[0].grad_gates[0].matrix,
+               1e-4);
+
+  GradientOfGate tmp2;
+  PopulateHessianFsimPhi("TheSymbol2", 0, 0, 1, 1.0, 1.0, 2.0, 1.0, &tmp2);
+
+  Matrix4Equal(tmp2.grad_gates[0].matrix, grad_gates[0].grad_gates[1].matrix,
+               1e-4);
+
+  metadata.clear();
+  meta.symbol_values.clear();
+  meta.placeholder_names.clear();
+  grad_gates.clear();
+  fuses.clear();
+
+  metadata.push_back(meta);
+
+  CreateHessianCircuit(circuit, metadata, &fuses, &grad_gates);
+  EXPECT_EQ(grad_gates.size(), 0);
+  EXPECT_EQ(fuses.size(), 1);
+}
+
+TEST(AdjHessianUtilTest, CreateHessianEmpty) {
+  QsimCircuit empty_circuit;
+  std::vector<GateMetaData> empty_metadata;
+  std::vector<std::vector<qsim::GateFused<QsimGate>>> fuses;
+  std::vector<GradientOfGate> grad_gates;
+
+  CreateHessianCircuit(empty_circuit, empty_metadata, &fuses, &grad_gates);
+
+  // Should create a single "empty fuse."
+  EXPECT_EQ(fuses.size(), 1);
+  EXPECT_EQ(fuses[0].size(), 0);
+
+  // No gradients.
+  EXPECT_EQ(grad_gates.size(), 0);
+}
+
+TEST(AdjHessianUtilTest, SingleEigenGrad) {
+  GradientOfGate grad;
+
+  PopulateHessianSingleEigen(&qsim::Cirq::YPowGate<float>::Create, "hello", 5,
+                             2, 0.125, 1.0, 0.0, &grad);
+
+  // Value verified from:
+  /*
+  (cirq.unitary(cirq.Y**(0.125 + 1e-2))
+    + cirq.unitary(cirq.Y**(0.125 - 1e-2))
+    - cirq.unitary(cirq.Y**(0.125))
+    - cirq.unitary(cirq.Y**(0.125))) * 1e4
+  array([[-4.55878779-1.88831173j,  1.88831173-4.55878779j],
+         [-1.88831173+4.55878779j, -4.55878779-1.88831173j]])
+  */
+  std::vector<float> expected{-4.558788,  -1.8883117, 1.8883117, -4.558788,
+                              -1.8883117, 4.558788,   -4.558788, -1.8883117};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hello");
+  Matrix2Equal(grad.grad_gates[0].matrix, expected, 5e-3);
+}
+
+TEST(AdjHessianUtilTest, TwoEigenGrad) {
+  GradientOfGate grad;
+
+  PopulateHessianTwoEigen(&qsim::Cirq::XXPowGate<float>::Create, "hi", 5, 2, 3,
+                          0.001, 1.0, 0.0, &grad);
+
+  // Value verified from:
+  /*
+  (cirq.unitary(cirq.XX**(0.001 + 1e-2))
+    + cirq.unitary(cirq.XX**(0.001 - 1e-2))
+    - cirq.unitary(cirq.XX**(0.001))
+    - cirq.unitary(cirq.XX**(0.001))) * 1e4
+    array([[-0.00049348-1.55031128e-06j,  0.        +0.00000000e+00j,
+             0.        +0.00000000e+00j,  0.00049348+1.55031128e-06j],
+           [ 0.        +0.00000000e+00j, -0.00049348-1.55031128e-06j,
+             0.00049348+1.55031128e-06j,  0.        +0.00000000e+00j],
+           [ 0.        +0.00000000e+00j,  0.00049348+1.55031128e-06j,
+            -0.00049348-1.55031128e-06j,  0.        +0.00000000e+00j],
+           [ 0.00049348+1.55031128e-06j,  0.        +0.00000000e+00j,
+             0.        +0.00000000e+00j, -0.00049348-1.55031128e-06j]])
+  */
+  std::vector<float> expected{
+      -4.934372, -0.01550184, 0.,        0.,          0.,        0.,
+      4.934372,  0.01550184,  0.,        0.,          -4.934372, -0.01550184,
+      4.934372,  0.01550184,  0.,        0.,          0.,        0.,
+      4.934372,  0.01550184,  -4.934372, -0.01550184, 0.,        0.,
+      4.934372,  0.01550184,  0.,        0.,          0.,        0.,
+      -4.934372, -0.01550184};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hi");
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 1e-3);
+}
+
+TEST(AdjHessianUtilTest, PhasedXPhasedExponent) {
+  GradientOfGate grad;
+
+  PopulateHessianPhasedXPhasedExponent("hello2", 5, 2, 0.001, 1.0, 1.0, 1.0,
+                                       0.0, &grad);
+  /* Value verified from:
+  (cirq.unitary(cirq.PhasedXPowGate(exponent=1.0,phase_exponent=0.001 + 1e-2))
+    + cirq.unitary(cirq.PhasedXPowGate(exponent=1.0,phase_exponent=0.001 -
+  1e-2))
+    - cirq.unitary(cirq.PhasedXPowGate(exponent=1.0,phase_exponent=0.001))
+    - cirq.unitary(cirq.PhasedXPowGate(exponent=1.0,phase_exponent=0.001))) *
+  1e4 array([[ 0.        +0.j        , -9.86874398+0.03100368j],
+           [-9.86874398-0.03100368j,  0.        +0.j        ]])
+  */
+  std::vector<float> expected{0.,        0.,          -9.868744, 0.03100368,
+                              -9.868744, -0.03100368, 0.,        0.};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hello2");
+  Matrix2Equal(grad.grad_gates[0].matrix, expected, 5e-4);
+}
+
+TEST(AdjHessianUtilTest, PhasedXExponent) {
+  GradientOfGate grad;
+
+  PopulateHessianPhasedXExponent("hello3", 5, 2, 10.123, 1.0, 0.789, 1.0, 0.0,
+                                 &grad);
+  /* Value verified from:
+  (cirq.unitary(cirq.PhasedXPowGate(exponent=0.789+1e-2,phase_exponent=10.123))
+    +
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789-1e-2,phase_exponent=10.123))
+    - cirq.unitary(cirq.PhasedXPowGate(exponent=0.789,phase_exponent=10.123))
+    - cirq.unitary(cirq.PhasedXPowGate(exponent=0.789,phase_exponent=10.123))) *
+  1e4 array([[ 3.88941758-3.03656025j, -2.45824276+4.2784705j ],
+           [-4.74702582+1.34685303j,  3.88941758-3.03656025j]])
+  */
+  std::vector<float> expected{3.8894176, -3.0365603, -2.4582427, 4.2784705,
+                              -4.747026, 1.346853,   3.8894176,  -3.0365603};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hello3");
+  Matrix2Equal(grad.grad_gates[0].matrix, expected, 1e-3);
+}
+
+TEST(AdjHessianUtilTest, CrossTermPhasedXPhasedExponentExponent) {
+  GradientOfGate grad;
+  // Symbol is always `kUsePrevTwoSymbols`, so this has no input argument for
+  // it.
+  PopulateCrossTermPhasedXPhasedExponentExponent(5, 2, 10.123, 1.0, 0.789, 1.0,
+                                                 0.0, &grad);
+  /* Value verified from:
+  (cirq.unitary(cirq.PhasedXPowGate(exponent=0.789+5e-3,phase_exponent=10.123+5e-3))
+    +
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789-5e-3,phase_exponent=10.123-5e-3))
+    -
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789+5e-3,phase_exponent=10.123-5e-3))
+    -
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789-5e-3,phase_exponent=10.123+5e-3)))
+  * 1e4 array([[
+  0.00000000e+00+5.55111512e-13j,  2.45824276e+00-4.27847050e+00j],
+           [-4.74702582e+00+1.34685303e+00j,  2.77555756e-13+5.55111512e-13j]])
+  */
+  std::vector<float> expected{0.0000000e+00,  5.5511151e-13,  2.4582427e+00,
+                              -4.2784705e+00, -4.7470260e+00, 1.3468530e+00,
+                              2.7755576e-13,  5.5511151e-13};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], kUsePrevTwoSymbols);
+  Matrix2Equal(grad.grad_gates[0].matrix, expected, 5e-3);
+}
+
+TEST(AdjHessianUtilTest, FSimThetaGrad) {
+  GradientOfGate grad;
+  PopulateHessianFsimTheta("hihi", 5, 2, 3, 0.5, 1.0, 1.2, 1.0, &grad);
+
+  /* Value verified from:
+  (cirq.unitary(cirq.FSimGate(theta=0.5 + 1e-2,phi=1.2))
+    + cirq.unitary(cirq.FSimGate(theta=0.5-1e-2,phi=1.2))
+    - cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2))
+    - cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2))) * 1e4
+    array([[ 0.        +0.j        ,  0.        +0.j        ,
+             0.        +0.j        ,  0.        +0.j        ],
+           [ 0.        +0.j        , -0.87757525+0.j        ,
+             0.        +0.47942154j,  0.        +0.j        ],
+           [ 0.        +0.j        ,  0.        +0.47942154j,
+             -0.87757525+0.j        ,  0.        +0.j        ],
+           [ 0.        +0.j        ,  0.        +0.j        ,
+            0.        +0.j        ,  0.        +0.j        ]])
+  */
+  std::vector<float> expected{
+      0., 0., 0.,         0.,         0.,         0.,         0., 0.,
+      0., 0., -0.8775753, 0.,         0.,         0.47942156, 0., 0.,
+      0., 0., 0.,         0.47942156, -0.8775753, 0.,         0., 0.,
+      0., 0., 0.,         0.,         0.,         0.,         0., 0.};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hihi");
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 1e-3);
+}
+
+TEST(AdjHessianUtilTest, FSimPhiGrad) {
+  GradientOfGate grad;
+  PopulateHessianFsimPhi("hihi2", 5, 2, 3, 0.5, 1.0, 1.2, 1.0, &grad);
+
+  /*
+  (cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2+1e-2))
+    + cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2-1e-2))
+    - cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2))
+    - cirq.unitary(cirq.FSimGate(theta=0.5,phi=1.2))) * 1e4
+  array([[ 0.        +0.j        ,  0.        +0.j        ,
+           0.        +0.j        ,  0.        +0.j        ],
+         [ 0.        +0.j        ,  0.        +0.j        ,
+           0.        +0.j        ,  0.        +0.j        ],
+         [ 0.        +0.j        ,  0.        +0.j        ,
+           0.        +0.j        ,  0.        +0.j        ],
+         [ 0.        +0.j        ,  0.        +0.j        ,
+           0.        +0.j        , -0.36235473+0.93203132j]])
+  */
+  std::vector<float> expected{0., 0., 0., 0., 0., 0., 0.,          0.,
+                              0., 0., 0., 0., 0., 0., 0.,          0.,
+                              0., 0., 0., 0., 0., 0., 0.,          0.,
+                              0., 0., 0., 0., 0., 0., -0.36235473, 0.93203133};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "hihi2");
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 5e-4);
+}
+
+TEST(AdjHessianUtilTest, CrossTermFsimThetaPhi) {
+  GradientOfGate grad;
+  // Symbol is always `kUsePrevTwoSymbols`, so this has no input argument for
+  // it.
+  PopulateCrossTermFsimThetaPhi(5, 2, 10.123, 1.0, 0.789, 1.0, 0.0, &grad);
+  /* Value verified from:
+  (cirq.unitary(cirq.PhasedXPowGate(exponent=0.789+5e-3,phase_exponent=10.123+5e-3))
+    +
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789-5e-3,phase_exponent=10.123-5e-3))
+    -
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789+5e-3,phase_exponent=10.123-5e-3))
+    -
+  cirq.unitary(cirq.PhasedXPowGate(exponent=0.789-5e-3,phase_exponent=10.123+5e-3)))
+  * 1e4 array([[0.00000000e+00+0.j, 0.00000000e+00+0.j, 0.00000000e+00+0.j,
+            0.00000000e+00+0.j],
+           [0.00000000e+00+0.j, 1.11022302e-12+0.j, 0.00000000e+00+0.j,
+            0.00000000e+00+0.j],
+           [0.00000000e+00+0.j, 0.00000000e+00+0.j, 1.11022302e-12+0.j,
+            0.00000000e+00+0.j],
+           [0.00000000e+00+0.j, 0.00000000e+00+0.j, 0.00000000e+00+0.j,
+            0.00000000e+00+0.j]])
+  */
+  std::vector<float> expected{0.,           0., 0., 0., 0., 0., 0., 0., 0., 0.,
+                              1.110223e-12, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+                              1.110223e-12, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+                              0.,           0.};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], kUsePrevTwoSymbols);
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 5e-3);
+}
+
+TEST(AdjHessianUtilTest, PhasedISwapPhasedExponent) {
+  GradientOfGate grad;
+
+  PopulateHessianPhasedISwapPhasedExponent("h", 5, 3, 2, 8.9, 1.0, -3.2, 1.0,
+                                           &grad);
+
+  /*
+  (cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9+1e-2))
+  + cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9-1e-2))
+  - cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9))
+  - cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9))) *
+  1e4
+
+  array([[  0.         +0.j        ,   0.         +0.j        ,
+            0.         +0.j        ,   0.         +0.j        ],
+         [  0.         +0.j        ,   0.         +0.j        ,
+            -22.06184686-30.36552715j,   0.         +0.j        ],
+         [  0.         +0.j        ,  22.06184686-30.36552715j,
+            0.         +0.j        ,   0.         +0.j        ],
+         [  0.         +0.j        ,   0.         +0.j        ,
+            0.         +0.j        ,   0.         +0.j        ]])
+  */
+  std::vector<float> expected{
+      0., 0.,         0.,         0., 0., 0., 0., 0.,        0.,         0., 0.,
+      0., -22.061848, -30.365528, 0., 0., 0., 0., 22.061848, -30.365528, 0., 0.,
+      0., 0.,         0.,         0., 0., 0., 0., 0.,        0.,         0.};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "h");
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 3e-2);
+}
+
+TEST(AdjHessianUtilTest, PhasedISwapExponent) {
+  GradientOfGate grad;
+
+  PopulateHessianPhasedISwapExponent("h2", 5, 3, 2, 8.9, 1.0, -3.2, 1.0, &grad);
+
+  /*
+  (cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2+1e-2,phase_exponent=8.9))
+  + cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2-1e-2,phase_exponent=8.9))
+  - cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9))
+  - cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9))) *
+  1e4 array([[ 0.        +0.j       ,  0.        +0.j       , 0.        +0.j ,
+  0.        +0.j       ], [ 0.        +0.j       , -0.76245319+0.j       ,
+          -1.37929079-1.8984309j,  0.        +0.j       ],
+         [ 0.        +0.j       ,  1.37929079-1.8984309j,
+          -0.76245319+0.j       ,  0.        +0.j       ],
+         [ 0.        +0.j       ,  0.        +0.j       ,
+           0.        +0.j       ,  0.        +0.j       ]])
+  */
+  std::vector<float> expected{
+      0., 0., 0.,         0.,        0.,         0.,        0., 0.,
+      0., 0., -0.7624532, 0.,        -1.3792908, -1.898431, 0., 0.,
+      0., 0., 1.3792908,  -1.898431, -0.7624532, 0.,        0., 0.,
+      0., 0., 0.,         0.,        0.,         0.,        0., 0.};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], "h2");
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 2e-3);
+}
+
+TEST(AdjHessianUtilTest, CrossTermPhasedISwapPhasedExponentExponent) {
+  GradientOfGate grad;
+  // Symbol is always `kUsePrevTwoSymbols`, so this has no input argument for
+  // it.
+  PopulateCrossTermPhasedISwapPhasedExponentExponent(5, 3, 2, 8.9, 1.0, -3.2,
+                                                     1.0, &grad);
+  /* Value verified from:
+  (cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2+5e-3,phase_exponent=8.9+5e-3))
+  +
+  cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2-5e-3,phase_exponent=8.9-5e-3))
+  -
+  cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2+5e-3,phase_exponent=8.9-5e-3))
+  -
+  cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2-5e-3,phase_exponent=8.9+5e-3)))
+  * 1e4 array([[ 0.00000000e+00+0.j        ,  0.00000000e+00+0.j        ,
+             0.00000000e+00+0.j        ,  0.00000000e+00+0.j        ],
+           [ 0.00000000e+00+0.j        , -5.55111512e-13+0.j        ,
+            -2.46696989e+00+1.79235854j,  0.00000000e+00+0.j        ],
+           [ 0.00000000e+00+0.j        ,  2.46696989e+00+1.79235854j,
+            -5.55111512e-13+0.j        ,  0.00000000e+00+0.j        ],
+           [ 0.00000000e+00+0.j        ,  0.00000000e+00+0.j        ,
+             0.00000000e+00+0.j        ,  0.00000000e+00+0.j        ]])
+  */
+  std::vector<float> expected{
+      0.0000000e+00,  0.0000000e+00, 0.0000000e+00,  0.0000000e+00,
+      0.0000000e+00,  0.0000000e+00, 0.0000000e+00,  0.0000000e+00,
+      0.0000000e+00,  0.0000000e+00, -5.5511151e-13, 0.0000000e+00,
+      -2.4669700e+00, 1.7923585e+00, 0.0000000e+00,  0.0000000e+00,
+      0.0000000e+00,  0.0000000e+00, 2.4669700e+00,  1.7923585e+00,
+      -5.5511151e-13, 0.0000000e+00, 0.0000000e+00,  0.0000000e+00,
+      0.0000000e+00,  0.0000000e+00, 0.0000000e+00,  0.0000000e+00,
+      0.0000000e+00,  0.0000000e+00, 0.0000000e+00,  0.0000000e+00};
+
+  EXPECT_EQ(grad.index, 5);
+  EXPECT_EQ(grad.params[0], kUsePrevTwoSymbols);
+  Matrix4Equal(grad.grad_gates[0].matrix, expected, 5e-3);
+}
+
+TEST(AdjHessianUtilTest, Matrix2Add) {
+  std::array<float, 8> u{1, 2, 3, 4, 5, 6, 7, 8};
+  std::array<float, 8> u2{0, -1, -2, -3, -4, -5, -6, -7};
+  Matrix2Add(u, u2);
+  for (int i = 0; i < 8; i++) {
+    EXPECT_EQ(u2[i], 1);
+    EXPECT_EQ(u[i], i + 1);
+  }
+}
+
+TEST(AdjHessianUtilTest, Matrix4Add) {
+  std::array<float, 32> u;
+  std::array<float, 32> u2;
+
+  for (int i = 0; i < 32; i++) {
+    u2[i] = -i;
+    u[i] = i + 1;
+  }
+
+  Matrix4Add(u, u2);
+  for (int i = 0; i < 32; i++) {
+    EXPECT_EQ(u2[i], 1);
+    EXPECT_EQ(u[i], i + 1);
+  }
+}
+
+}  // namespace
+}  // namespace tfq
diff --git a/tensorflow_quantum/core/src/adj_util_test.cc b/tensorflow_quantum/core/src/adj_util_test.cc
index 526efc175..a79bca1b3 100644
--- a/tensorflow_quantum/core/src/adj_util_test.cc
+++ b/tensorflow_quantum/core/src/adj_util_test.cc
@@ -507,7 +507,8 @@ TEST(AdjUtilTest, PhasedISwapPhasedExponent) {
 
   /*
   (cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9+1e-4))
-  - cirq.unitary(cirq.PhasedISwapPowGate(exponent=3.2,phase_exponent=8.9-1e-4)))
+  -
+  cirq.unitary(cirq.PhasedISwapPowGate(exponent=-3.2,phase_exponent=8.9-1e-4)))
     / 2e-4
   array([[ 0.        +0.j        ,  0.        +0.j        ,
            0.        +0.j        ,  0.        +0.j        ],