diff --git a/docs/sphinx/examples/python/tutorials/hadamard_test.ipynb b/docs/sphinx/examples/python/tutorials/hadamard_test.ipynb
new file mode 100644
index 0000000000..2c811ad7a9
--- /dev/null
+++ b/docs/sphinx/examples/python/tutorials/hadamard_test.ipynb
@@ -0,0 +1,344 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Hadamard Test\n",
+    "\n",
+    "Consider the observable $O$ and two generic quantum states $\\ket{\\psi}$ and $\\ket{\\phi}$. We want to calculate the quantity\n",
+    "$$\n",
+    "\\bra{\\psi} O \\ket{\\phi}.\n",
+    "$$\n",
+    "where $O$ is a Pauli operator.\n",
+    "\n",
+    "First of all we shall prepare the states $\\ket{\\psi}$ and $\\ket{\\phi}$ using a quantum circuit for each of them. So we  have\n",
+    "$$\n",
+    "\\ket{\\psi} = U_{\\psi}\\ket{0} \\qquad \\ket{\\phi} = U_{\\phi}\\ket{0}\n",
+    "$$\n",
+    "\n",
+    "Let's define an observable we want to use:\n",
+    "$$\n",
+    "O = X_1X_2\n",
+    "$$\n",
+    "\n",
+    "Now we can evaluate the matrix element using the following fact:\n",
+    "$$\n",
+    "\\bra{\\psi}O\\ket{\\phi} = \\bra{0}U_\\psi^\\dagger O U_\\phi\\ket{0}\n",
+    "$$\n",
+    "This is just an expectation value which can be solved with a simple Hadamard test. The probability to measure $0$ or $1$ in the ancilla qubit is\n",
+    "\n",
+    "$$\n",
+    "P(0) = \\frac{1}{2} \\left[ I + Re \\bra{\\psi} O \\ket{\\phi} \\right]\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "P(1) = \\frac{1}{2} \\left[ I - Re \\bra{\\psi} O \\ket{\\phi} \\right]\n",
+    "$$\n",
+    "\n",
+    "The difference between the probability of $0$ and $1$ gives \n",
+    "\n",
+    "$$\n",
+    "P(0)-P(1) = Re \\bra{\\psi} O \\ket{\\phi}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### A- Numerical result as a reference: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Psi state:  (0.707107,0)\n",
+      "(0,0)\n",
+      "(0.707107,0)\n",
+      "(0,0)\n",
+      "\n",
+      "Phi state:  (0,0)\n",
+      "(1,0)\n",
+      "(0,0)\n",
+      "(0,0)\n",
+      "\n",
+      "hamiltonian:  [[0.+0.j 0.+0.j 0.+0.j 1.+0.j]\n",
+      " [0.+0.j 0.+0.j 1.+0.j 0.+0.j]\n",
+      " [0.+0.j 1.+0.j 0.+0.j 0.+0.j]\n",
+      " [1.+0.j 0.+0.j 0.+0.j 0.+0.j]] \n",
+      "\n",
+      "Numerical expectation value:  (0.7071067690849304+0j)\n"
+     ]
+    }
+   ],
+   "source": [
+    "import cudaq\n",
+    "import numpy as np\n",
+    "from functools import reduce\n",
+    "\n",
+    "cudaq.set_target('nvidia')\n",
+    "\n",
+    "qubit_num = 2\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def psi(num:int):\n",
+    "    q = cudaq.qvector(num)\n",
+    "    h(q[1])\n",
+    "    \n",
+    "@cudaq.kernel\n",
+    "def phi(n:int):\n",
+    "    q = cudaq.qvector(n)\n",
+    "    x(q[0])\n",
+    "\n",
+    "psi_state = cudaq.get_state(psi, qubit_num)\n",
+    "print('Psi state: ', psi_state)\n",
+    "\n",
+    "phi_state=cudaq.get_state(phi, qubit_num)\n",
+    "print('Phi state: ', phi_state)\n",
+    "\n",
+    "ham=cudaq.spin.x(0) * cudaq.spin.x(1)\n",
+    "ham_matrix = ham.to_matrix()\n",
+    "print('hamiltonian: ', np.array(ham_matrix), '\\n')\n",
+    "\n",
+    "exp_val=reduce(np.dot,(np.array(psi_state).conj().T, ham_matrix, phi_state))\n",
+    "\n",
+    "print('Numerical expectation value: ', exp_val) "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### B- Using ``sample`` algorithmic primitive to sample the ancilla qubit and compute the expectation value."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{ 0:85356 1:14644 }\n",
+      "\n",
+      "Observable QC:  0.70712 + - 0.0015811092713661505\n",
+      "Numerical result 0.7071067690849304\n"
+     ]
+    }
+   ],
+   "source": [
+    "import cudaq\n",
+    "\n",
+    "cudaq.set_target('nvidia')\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def U_psi(q:cudaq.qview):\n",
+    "    h(q[1])\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def U_phi(q:cudaq.qview):\n",
+    "    x(q[0])\n",
+    "\n",
+    "@cudaq.kernel  \n",
+    "def ham_cir(q:cudaq.qview):\n",
+    "    x(q[0])\n",
+    "    x(q[1])\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def kernel(n:int):\n",
+    "    ancilla=cudaq.qubit()\n",
+    "    q = cudaq.qvector(n)\n",
+    "    h(ancilla)\n",
+    "    cudaq.control(U_phi,ancilla,q)\n",
+    "    cudaq.control(ham_cir,ancilla,q)\n",
+    "    cudaq.control(U_psi,ancilla,q)\n",
+    "    \n",
+    "    h(ancilla)\n",
+    "    \n",
+    "    mz(ancilla)\n",
+    "\n",
+    "shots = 100000    \n",
+    "qubit_num=2\n",
+    "count = cudaq.sample(kernel, qubit_num, shots_count = shots)    \n",
+    "print(count)\n",
+    "\n",
+    "mean_val = (count['0']-count['1']) / shots\n",
+    "error = np.sqrt(2* count['0'] * count['1'] / shots) / shots\n",
+    "print('Observable QC: ', mean_val,'+ -', error)\n",
+    "print('Numerical result', np.real(exp_val))\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### C- Use multi-GPUs to compute the matrix elements"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Number of QPUs: 5\n",
+      "0\n",
+      "{ 0:63807 1:36193 }\n",
+      "\n",
+      "Observable QC:  0.27614 + - 0.0021491238917289066\n",
+      "1\n",
+      "{ 0:49929 1:50071 }\n",
+      "\n",
+      "Observable QC:  -0.00142 + - 0.0022360657230949183\n",
+      "2\n",
+      "{ 0:50041 1:49959 }\n",
+      "\n",
+      "Observable QC:  0.00082 + - 0.0022360672257336093\n",
+      "3\n",
+      "{ 0:50276 1:49724 }\n",
+      "\n",
+      "Observable QC:  0.00552 + - 0.0022360339102974265\n"
+     ]
+    }
+   ],
+   "source": [
+    "import cudaq\n",
+    "\n",
+    "cudaq.set_target(\"nvidia-mqpu\")\n",
+    "\n",
+    "target = cudaq.get_target()\n",
+    "qpu_count = target.num_qpus()\n",
+    "print(\"Number of QPUs:\", qpu_count)\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def U_psi(q:cudaq.qview, theta:float):\n",
+    "    ry(theta, q[1])\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def U_phi(q:cudaq.qview, theta: float):\n",
+    "    rx(theta, q[0])\n",
+    "\n",
+    "@cudaq.kernel  \n",
+    "def ham_cir(q:cudaq.qview):\n",
+    "    x(q[0])\n",
+    "    x(q[1])\n",
+    "\n",
+    "@cudaq.kernel\n",
+    "def kernel(n:int, angle:float, theta:float):\n",
+    "    ancilla = cudaq.qubit()\n",
+    "    q = cudaq.qvector(n)\n",
+    "    h(ancilla)\n",
+    "    cudaq.control(U_phi, ancilla, q, theta)\n",
+    "    cudaq.control(ham_cir, ancilla, q)\n",
+    "    cudaq.control(U_psi, ancilla, q, angle)\n",
+    "    \n",
+    "    h(ancilla)\n",
+    "    \n",
+    "    mz(ancilla)\n",
+    "    \n",
+    "shots = 100000  \n",
+    "angle = [0.0, 1.5,3.14,0.7]\n",
+    "theta = [0.6, 1.2 ,2.2 ,3.0]\n",
+    "qubit_num = 2\n",
+    "\n",
+    "result = []\n",
+    "for i in range(4):  \n",
+    "    count = cudaq.sample_async(kernel, qubit_num, angle[i], theta[i], shots_count = shots, qpu_id = i)  \n",
+    "    result.append(count)  \n",
+    "\n",
+    "mean_val = np.zeros(len(angle))\n",
+    "i = 0\n",
+    "for count in result:\n",
+    "    print(i)\n",
+    "    i_result = count.get()\n",
+    "    print(i_result)\n",
+    "    mean_val[i] = (i_result['0'] - i_result['1']) / shots\n",
+    "    error = np.sqrt(2 * i_result['0'] * i_result['1'] / shots) / shots\n",
+    "    print('Observable QC: ',  mean_val[i],'+ -', error)\n",
+    "    i += 1\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Diagonalize the matrix using for example Numpy or CuPy. In this example, since we are having 2x2 matrix, we use numpy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 0.27614 -0.00142]\n",
+      " [ 0.00082  0.00552]]\n",
+      "Eigen values: \n",
+      "[0.00551752 0.27614248]\n",
+      "Eigenvector: \n",
+      "[[ 0.00303004 -0.99999541]\n",
+      " [-0.99999541 -0.00303004]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "my_mat = np.zeros((2,2),dtype=float)\n",
+    "m = 0\n",
+    "for k in range(2):\n",
+    "    for j in range(2):\n",
+    "        my_mat[k,j] = mean_val[m]\n",
+    "        m += 1   \n",
+    "\n",
+    "print(my_mat)\n",
+    "\n",
+    "E,V = np.linalg.eigh(my_mat)\n",
+    "\n",
+    "print('Eigen values: ')\n",
+    "print(E)\n",
+    "\n",
+    "print('Eigenvector: ')\n",
+    "print(V)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/sphinx/using/tutorials.rst b/docs/sphinx/using/tutorials.rst
index 14fcf89279..6a44187e8f 100644
--- a/docs/sphinx/using/tutorials.rst
+++ b/docs/sphinx/using/tutorials.rst
@@ -11,6 +11,7 @@ Tutorials that give an in depth view of CUDA-Q and its applications in Python.
     /examples/python/tutorials/cost_minimization.ipynb
     /examples/python/tutorials/vqe.ipynb
     /examples/python/tutorials/qaoa.ipynb
+    /examples/python/tutorials/hadamard_test.ipynb
     /examples/python/tutorials/hybrid_qnns.ipynb
     /examples/python/tutorials/noisy_simulations.ipynb
     /examples/python/tutorials/readout_error_mitigation.ipynb