diff --git a/.github/workflows/config/spelling_allowlist.txt b/.github/workflows/config/spelling_allowlist.txt
index 3fe9cc4537..8d88037d31 100644
--- a/.github/workflows/config/spelling_allowlist.txt
+++ b/.github/workflows/config/spelling_allowlist.txt
@@ -258,6 +258,8 @@ qubit
 qubits
 qudit
 qudits
+qumode
+qumodes
 reStructuredText
 reconfigurable
 runtime
diff --git a/docs/sphinx/api/default_ops.rst b/docs/sphinx/api/default_ops.rst
index 62d239b421..ac1d70b86c 100644
--- a/docs/sphinx/api/default_ops.rst
+++ b/docs/sphinx/api/default_ops.rst
@@ -638,31 +638,67 @@ operations, each operating on 2 qubits.
 Photonic Operations on Qudits
 =============================
 
-These operations are valid only on the `orca-photonics` target which does not support the quantum operations above.
+These operations are valid only on the `orca-photonics` target which does not support
+the quantum operations above.
 
-:code:`plus`
+:code:`create`
 ---------------------
 
-This is a place-holder, to be updated later.
+This operation increments the number of photons in a qumode up to a maximum value
+defined by the qudit level that represents the qumode. If it is applied to a qumode
+where the number of photons is already at the maximum value, the operation has no
+effect.
+
+:math:`U|0\rangle → |1\rangle, U|1\rangle → |2\rangle, U|2\rangle → |3\rangle, \cdots, U|d\rangle → |d\rangle`
+where :math:`d` is the qudit level.
 
 .. tab:: Python
 
     .. code-block:: python
 
         q = qudit(3)
-        plus(q)
+        create(q)
 
 .. tab:: C++
 
     .. code-block:: cpp
 
         cudaq::qvector<3> q(1);
-        plus(q[0]);
+        create(q[0]);
+
+:code:`annihilate`
+---------------------
+
+This operation reduces the number of photons in a qumode up to a minimum value of
+0 representing the vacuum state. If it is applied to a qumode where the number of
+photons is already at the minimum value 0, the operation has no effect.
+
+:math:`U|0\rangle → |0\rangle, U|1\rangle → |0\rangle, U|2\rangle → |1\rangle, \cdots, U|d\rangle → |d-1\rangle`
+where :math:`d` is the qudit level.
+
+.. tab:: Python
+
+    .. code-block:: python
+
+        q = qudit(3)
+        annihilate(q)
+
+.. tab:: C++
+
+    .. code-block:: cpp
+
+        cudaq::qvector<3> q(1);
+        annihilate(q[0]);
 
 :code:`phase_shift`
 ---------------------
 
-This is a place-holder, to be updated later.
+A phase shifter adds a phase :math:`\phi` on a qumode. For the annihilation (:math:`a_1`)
+and creation operators (:math:`a_1^\dagger`) of a qumode, the phase shift operator
+is defined  by
+
+.. math::
+    P(\phi) = \exp\left(i \phi a_1^\dagger a_1  \right)
 
 .. tab:: Python
 
@@ -681,7 +717,13 @@ This is a place-holder, to be updated later.
 :code:`beam_splitter`
 ---------------------
 
-This is a place-holder, to be updated later.
+Beam splitters act on two qumodes together and it is parameterized by a single angle 
+:math:`\theta`, relating to reflectivity.
+For the annihilation (:math:`a_1` and :math:`a_2`) and creation operators (:math:`a_1^\dagger`
+and :math:`a_2^\dagger`) of two qumodes, the beam splitter operator is defined by
+
+.. math::
+    B(\theta) = \exp\left[i \theta (a_1^\dagger a_2 + a_1 a_2^\dagger) \right]
 
 .. tab:: Python
 
@@ -700,19 +742,19 @@ This is a place-holder, to be updated later.
 :code:`mz`
 ---------------------
 
-This operation returns the measurement results of the input qudit(s).
+This operation returns the measurement results of the input qumode(s).
 
 .. tab:: Python
 
     .. code-block:: python
 
-        qutrits = [qudit(3) for _ in range(2)]
-        mz(qutrits)
+        qumodes = [qudit(3) for _ in range(2)]
+        mz(qumodes)
 
 
 .. tab:: C++
 
     .. code-block:: cpp
 
-        cudaq::qvector<3> qutrits(2);
-        mz(qutrits);
+        cudaq::qvector<3> qumodes(2);
+        mz(qumodes);
diff --git a/docs/sphinx/targets/cpp/orca.cpp b/docs/sphinx/targets/cpp/orca.cpp
index 289145ee34..78f8a00a57 100644
--- a/docs/sphinx/targets/cpp/orca.cpp
+++ b/docs/sphinx/targets/cpp/orca.cpp
@@ -5,24 +5,13 @@
 // To use the ORCA Computing target you will need to set the ORCA_ACCESS_URL
 // environment variable or pass the URL to the `--orca-url` flag.
 
-#include "cudaq/orca.h"
-#include "cudaq.h"
 #include <chrono>
+#include <cudaq.h>
+#include <cudaq/orca.h>
 #include <fstream>
 #include <iostream>
 #include <thread>
 
-// define helper function to generate linear spaced vectors
-template <typename T>
-void linear_spaced_vector(std::vector<T> &xs, T min, T max, std::size_t N) {
-  T h = (max - min) / static_cast<T>(N - 1);
-  typename std::vector<T>::iterator x;
-  T val;
-  for (x = xs.begin(), val = min; x != xs.end(); ++x, val += h) {
-    *x = val;
-  }
-}
-
 int main() {
   using namespace std::this_thread;     // sleep_for, sleep_until
   using namespace std::chrono_literals; // `ns`, `us`, `ms`, `s`, `h`, etc.
@@ -39,10 +28,10 @@ int main() {
 
   // half of 8 time bins is filled with a single photon and the other half is
   // filled with the vacuum state (empty)
-  std::vector<std::size_t> input_state{1, 0, 1, 0, 1, 0, 1, 0};
+  std::vector<std::size_t> input_state = {1, 0, 1, 0, 1, 0, 1, 0};
 
   // The time bin interferometer in this example has two loops, each of length 1
-  std::vector<std::size_t> loop_lengths{1, 1};
+  std::vector<std::size_t> loop_lengths = {1, 1};
 
   // helper variables to calculate the number of beam splitters and phase
   // shifters needed in the TBI
@@ -53,18 +42,18 @@ int main() {
   const std::size_t n_beam_splitters = n_loops * n_modes - sum_loop_lengths;
 
   // beam splitter angles (created as a linear spaced vector of angles)
-  std::vector<double> bs_angles(n_beam_splitters);
-  linear_spaced_vector(bs_angles, M_PI / 8, M_PI / 3, n_beam_splitters);
+  std::vector<double> bs_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 6, n_beam_splitters);
 
   // Optionally, we can also specify the phase shifter angles (created as a
   // linear spaced vector of angles), if the system includes phase shifters
   // ```
-  // std::vector<double> ps_angles(n_beam_splitters);
-  // linear_spaced_vector(ps_angles, M_PI / 6, M_PI / 3, n_beam_splitters);
+  // std::vector<double> ps_angles = cudaq::linspace(M_PI / 3, M_PI / 5,
+  // n_beam_splitters);
   // ```
 
   // we can also set number of requested samples
-  int n_samples{10000};
+  int n_samples = 10000;
 
   // Submit to ORCA synchronously (e.g., wait for the job result to be
   // returned before proceeding with the rest of the execution).
diff --git a/docs/sphinx/targets/cpp/orca_mqpu.cpp b/docs/sphinx/targets/cpp/orca_mqpu.cpp
index 496c4454f8..96b3f68b1b 100644
--- a/docs/sphinx/targets/cpp/orca_mqpu.cpp
+++ b/docs/sphinx/targets/cpp/orca_mqpu.cpp
@@ -5,23 +5,11 @@
 // ```
 // See accompanying example `orca.cpp` for detailed explanation.
 
-#include "cudaq.h"
-#include "cudaq/orca.h"
-
+#include <cudaq.h>
+#include <cudaq/orca.h>
 #include <fstream>
 #include <iostream>
 
-// define helper function to generate linear spaced vectors
-template <typename T>
-void linear_spaced_vector(std::vector<T> &xs, T min, T max, std::size_t N) {
-  T h = (max - min) / static_cast<T>(N - 1);
-  typename std::vector<T>::iterator x;
-  T val;
-  for (x = xs.begin(), val = min; x != xs.end(); ++x, val += h) {
-    *x = val;
-  }
-}
-
 int main() {
 
   auto &platform = cudaq::get_platform();
@@ -29,16 +17,16 @@ int main() {
   printf("Number of QPUs: %zu\n", num_qpus);
 
   // A time-bin boson sampling experiment
-  std::vector<std::size_t> input_state{1, 0, 1, 0, 1, 0, 1, 0};
-  std::vector<std::size_t> loop_lengths{1, 1};
+  std::vector<std::size_t> input_state = {1, 0, 1, 0, 1, 0, 1, 0};
+  std::vector<std::size_t> loop_lengths = {1, 1};
   std::size_t sum_loop_lengths{std::accumulate(
       loop_lengths.begin(), loop_lengths.end(), static_cast<std::size_t>(0))};
   const std::size_t n_loops = loop_lengths.size();
   const std::size_t n_modes = input_state.size();
   const std::size_t n_beam_splitters = n_loops * n_modes - sum_loop_lengths;
-  std::vector<double> bs_angles(n_beam_splitters);
-  linear_spaced_vector(bs_angles, M_PI / 8, M_PI / 3, n_beam_splitters);
-  int n_samples{10000};
+  std::vector<double> bs_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 6, n_beam_splitters);
+  int n_samples = 10000;
 
   std::cout << "Submitting to ORCA Server asynchronously" << std::endl;
   std::vector<cudaq::async_sample_result> countFutures;
diff --git a/docs/sphinx/targets/cpp/photonics.cpp b/docs/sphinx/targets/cpp/photonics.cpp
index df81ad3a50..bc69ee66be 100644
--- a/docs/sphinx/targets/cpp/photonics.cpp
+++ b/docs/sphinx/targets/cpp/photonics.cpp
@@ -9,20 +9,18 @@
 
 struct photonicsKernel {
   void operator()() __qpu__ {
-    cudaq::qvector<3> qutrits(2);
-    plus(qutrits[0]);
-    plus(qutrits[1]);
-    plus(qutrits[1]);
-    mz(qutrits);
+    cudaq::qvector<3> qumodes(2);
+    create(qumodes[0]);
+    create(qumodes[1]);
+    create(qumodes[1]);
+    mz(qumodes);
   }
 };
 
 int main() {
 
   auto counts = cudaq::sample(photonicsKernel{});
-  for (auto &[k, v] : counts) {
-    printf("Result : Count = %s : %lu\n", k.c_str(), v);
-  }
+  counts.dump();
 
   auto state = cudaq::get_state(photonicsKernel{});
   state.dump();
diff --git a/docs/sphinx/targets/cpp/photonics_tbi.cpp b/docs/sphinx/targets/cpp/photonics_tbi.cpp
deleted file mode 100644
index d9c5701b06..0000000000
--- a/docs/sphinx/targets/cpp/photonics_tbi.cpp
+++ /dev/null
@@ -1,84 +0,0 @@
-// Compile and run with:
-// ```
-// nvq++ --target orca-photonics photonics_tbi.cpp -o tbi.x
-// ./tbi.x
-// ```
-
-#include "cudaq.h"
-#include "cudaq/photonics.h"
-
-#include <iostream>
-
-// Global variables
-static const std::size_t one{1};
-
-static constexpr std::size_t n_modes{4};
-static constexpr std::array<std::size_t, n_modes> input_state{2, 1, 3, 1};
-
-static constexpr std::size_t d{
-    std::accumulate(input_state.begin(), input_state.end(), one)};
-
-struct TBI {
-  auto operator()(std::vector<double> const &bs_angles,
-                  std::vector<double> const &ps_angles,
-                  std::vector<std::size_t> const &input_state,
-                  std::vector<std::size_t> const &loop_lengths) __qpu__ {
-    auto n_modes = ::n_modes;
-    const auto d = ::d;
-
-    cudaq::qvector<d> quds(n_modes); // |00...00> d-dimensions
-    for (std::size_t i = 0; i < n_modes; i++) {
-      for (std::size_t j = 0; j < input_state[i]; j++) {
-        plus(quds[i]); // setting to |input_state>
-      }
-    }
-
-    std::size_t c = 0;
-    for (std::size_t ll : loop_lengths) {
-      for (std::size_t i = 0; i < (n_modes - ll); i++) {
-        beam_splitter(quds[i], quds[i + ll], bs_angles[c]);
-        phase_shift(quds[i], ps_angles[c]);
-        c++;
-      }
-    }
-    mz(quds);
-  }
-};
-
-template <typename T>
-void LinearSpacedArray(std::vector<T> &xs, T min, T max, std::size_t N) {
-  T h = (max - min) / static_cast<T>(N - 1);
-  typename std::vector<T>::iterator x;
-  T val;
-  for (x = xs.begin(), val = min; x != xs.end(); ++x, val += h) {
-    *x = val;
-  }
-}
-
-int main() {
-  std::size_t n_loops{2};
-  std::vector<std::size_t> loop_lengths{1, 2};
-  std::vector<std::size_t> input_state(std::begin(::input_state),
-                                       std::end(::input_state));
-
-  const std::size_t zero{0};
-  std::size_t sum_loop_lenghts{
-      std::accumulate(loop_lengths.begin(), loop_lengths.end(), zero)};
-
-  std::size_t n_beam_splitters = n_loops * ::n_modes - sum_loop_lenghts;
-
-  std::vector<double> bs_angles(n_beam_splitters);
-  std::vector<double> ps_angles(n_beam_splitters);
-
-  LinearSpacedArray(bs_angles, M_PI / 3, M_PI / 6, n_beam_splitters);
-  LinearSpacedArray(ps_angles, M_PI / 3, M_PI / 5, n_beam_splitters);
-
-  auto counts = cudaq::sample(1000000, TBI{}, bs_angles, ps_angles, input_state,
-                              loop_lengths);
-
-  for (auto &[k, v] : counts) {
-    std::cout << k << ":" << v << " ";
-  }
-  std::cout << std::endl;
-  return 0;
-}
\ No newline at end of file
diff --git a/docs/sphinx/targets/cpp/photonics_tbi_get_state.cpp b/docs/sphinx/targets/cpp/photonics_tbi_get_state.cpp
new file mode 100644
index 0000000000..3b27be44f5
--- /dev/null
+++ b/docs/sphinx/targets/cpp/photonics_tbi_get_state.cpp
@@ -0,0 +1,69 @@
+// Compile and run with:
+// ```
+// nvq++ --target orca-photonics photonics_tbi_get_state.cpp && ./a.out
+// ```
+
+#include <cudaq.h>
+#include <cudaq/photonics.h>
+#include <iostream>
+
+// Global variables
+static constexpr std::size_t one = 1;
+
+static constexpr std::size_t n_modes = 8;
+static constexpr std::array<std::size_t, n_modes> input_state = {1, 0, 1, 0,
+                                                                 1, 0, 1, 0};
+
+static constexpr std::size_t d =
+    std::accumulate(input_state.begin(), input_state.end(), one);
+
+struct TBI {
+  auto operator()(std::vector<double> const &bs_angles,
+                  std::vector<double> const &ps_angles,
+                  std::vector<std::size_t> const &input_state,
+                  std::vector<std::size_t> const &loop_lengths) __qpu__ {
+    auto n_modes = ::n_modes;
+    const auto d = ::d;
+
+    cudaq::qvector<d> qumodes(n_modes); // |00...00> d-dimensions
+    for (std::size_t i = 0; i < n_modes; i++) {
+      for (std::size_t j = 0; j < input_state[i]; j++) {
+        create(qumodes[i]); // setting to |input_state>
+      }
+    }
+
+    std::size_t c = 0;
+    for (std::size_t ll : loop_lengths) {
+      for (std::size_t i = 0; i < (n_modes - ll); i++) {
+        beam_splitter(qumodes[i], qumodes[i + ll], bs_angles[c]);
+        phase_shift(qumodes[i], ps_angles[c]);
+        c++;
+      }
+    }
+  }
+};
+
+int main() {
+  std::size_t n_loops = 2;
+  std::vector<std::size_t> loop_lengths = {1, 1};
+  std::vector<std::size_t> input_state(std::begin(::input_state),
+                                       std::end(::input_state));
+
+  const std::size_t zero = 0;
+  std::size_t sum_loop_lenghts{
+      std::accumulate(loop_lengths.begin(), loop_lengths.end(), zero)};
+
+  std::size_t n_beam_splitters = n_loops * ::n_modes - sum_loop_lenghts;
+
+  std::vector<double> bs_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 6, n_beam_splitters);
+  std::vector<double> ps_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 5, n_beam_splitters);
+
+  auto state =
+      cudaq::get_state(TBI{}, bs_angles, ps_angles, input_state, loop_lengths);
+
+  state.dump();
+
+  return 0;
+}
\ No newline at end of file
diff --git a/docs/sphinx/targets/cpp/photonics_tbi_sample.cpp b/docs/sphinx/targets/cpp/photonics_tbi_sample.cpp
new file mode 100644
index 0000000000..90da940ec2
--- /dev/null
+++ b/docs/sphinx/targets/cpp/photonics_tbi_sample.cpp
@@ -0,0 +1,70 @@
+// Compile and run with:
+// ```
+// nvq++ --target orca-photonics photonics_tbi_sample.cpp && ./a.out
+// ```
+
+#include <cudaq.h>
+#include <cudaq/photonics.h>
+#include <iostream>
+
+// Global variables
+static constexpr std::size_t one = 1;
+
+static constexpr std::size_t n_modes = 8;
+static constexpr std::array<std::size_t, n_modes> input_state = {1, 0, 1, 0,
+                                                                 1, 0, 1, 0};
+
+static constexpr std::size_t d =
+    std::accumulate(input_state.begin(), input_state.end(), one);
+
+struct TBI {
+  auto operator()(std::vector<double> const &bs_angles,
+                  std::vector<double> const &ps_angles,
+                  std::vector<std::size_t> const &input_state,
+                  std::vector<std::size_t> const &loop_lengths) __qpu__ {
+    auto n_modes = ::n_modes;
+    const auto d = ::d;
+
+    cudaq::qvector<d> qumodes(n_modes); // |00...00> d-dimensions
+    for (std::size_t i = 0; i < n_modes; i++) {
+      for (std::size_t j = 0; j < input_state[i]; j++) {
+        create(qumodes[i]); // setting to |input_state>
+      }
+    }
+
+    std::size_t c = 0;
+    for (std::size_t ll : loop_lengths) {
+      for (std::size_t i = 0; i < (n_modes - ll); i++) {
+        beam_splitter(qumodes[i], qumodes[i + ll], bs_angles[c]);
+        phase_shift(qumodes[i], ps_angles[c]);
+        c++;
+      }
+    }
+    mz(qumodes);
+  }
+};
+
+int main() {
+  std::size_t n_loops = 2;
+  std::vector<std::size_t> loop_lengths = {1, 1};
+  std::vector<std::size_t> input_state(std::begin(::input_state),
+                                       std::end(::input_state));
+
+  const std::size_t zero = 0;
+  std::size_t sum_loop_lenghts{
+      std::accumulate(loop_lengths.begin(), loop_lengths.end(), zero)};
+
+  std::size_t n_beam_splitters = n_loops * ::n_modes - sum_loop_lenghts;
+
+  std::vector<double> bs_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 6, n_beam_splitters);
+  std::vector<double> ps_angles =
+      cudaq::linspace(M_PI / 3, M_PI / 5, n_beam_splitters);
+
+  auto counts = cudaq::sample(1000000, TBI{}, bs_angles, ps_angles, input_state,
+                              loop_lengths);
+
+  counts.dump();
+
+  return 0;
+}
\ No newline at end of file
diff --git a/docs/sphinx/targets/python/orca.py b/docs/sphinx/targets/python/orca.py
index 53cadb09d5..79f5cfd770 100644
--- a/docs/sphinx/targets/python/orca.py
+++ b/docs/sphinx/targets/python/orca.py
@@ -35,12 +35,12 @@
 n_beam_splitters = len(loop_lengths) * len(input_state) - sum(loop_lengths)
 
 # beam splitter angles
-bs_angles = np.linspace(np.pi / 8, np.pi / 3, n_beam_splitters)
+bs_angles = np.linspace(np.pi / 3, np.pi / 6, n_beam_splitters)
 
 # Optionally, we can also specify the phase shifter angles, if the system
 # includes phase shifters
 # ```
-# ps_angles = np.linspace(np.pi / 6, np.pi / 3, n_beam_splitters)
+# ps_angles = np.linspace(np.pi / 3, np.pi / 5, n_beam_splitters)
 # ```
 
 # we can also set number of requested samples
diff --git a/docs/sphinx/targets/python/orca_mqpu.py b/docs/sphinx/targets/python/orca_mqpu.py
index 869ae656bf..74fe7c1b11 100644
--- a/docs/sphinx/targets/python/orca_mqpu.py
+++ b/docs/sphinx/targets/python/orca_mqpu.py
@@ -18,7 +18,7 @@
 input_state = [1, 0, 1, 0, 1, 0, 1, 0]
 loop_lengths = [1, 1]
 n_beam_splitters = len(loop_lengths) * len(input_state) - sum(loop_lengths)
-bs_angles = np.linspace(np.pi / 8, np.pi / 3, n_beam_splitters)
+bs_angles = np.linspace(np.pi / 3, np.pi / 6, n_beam_splitters)
 n_samples = 10000
 
 count_futures = []
diff --git a/docs/sphinx/targets/python/photonics.py b/docs/sphinx/targets/python/photonics.py
index b11857b451..ecf1696f3f 100644
--- a/docs/sphinx/targets/python/photonics.py
+++ b/docs/sphinx/targets/python/photonics.py
@@ -6,9 +6,9 @@
 @cudaq.kernel
 def photonicsKernel():
     qutrits = [qudit(3) for _ in range(2)]
-    plus(qutrits[0])
-    plus(qutrits[1])
-    plus(qutrits[1])
+    create(qutrits[0])
+    create(qutrits[1])
+    create(qutrits[1])
     mz(qutrits)
 
 
diff --git a/docs/sphinx/targets/python/photonics_tbi_get_state.py b/docs/sphinx/targets/python/photonics_tbi_get_state.py
new file mode 100644
index 0000000000..6b341e58ce
--- /dev/null
+++ b/docs/sphinx/targets/python/photonics_tbi_get_state.py
@@ -0,0 +1,38 @@
+import cudaq
+import numpy as np
+
+cudaq.set_target("photonics")
+
+
+@cudaq.kernel
+def TBI(
+    bs_angles: list[float],
+    ps_angles: list[float],
+    input_state: list[int],
+    loop_lengths: list[int],
+):
+    n_modes = len(input_state)
+    level = sum(input_state) + 1  # qudit level
+
+    qumodes = [qudit(level) for _ in range(n_modes)]
+
+    for i in range(n_modes):
+        for _ in range(input_state[i]):
+            create(qumodes[i])
+
+    counter = 0
+    for j in loop_lengths:
+        for i in range(n_modes - j):
+            beam_splitter(qumodes[i], qumodes[i + j], bs_angles[counter])
+            phase_shift(qumodes[i], ps_angles[counter])
+            counter += 1
+
+
+input_state = [1, 0, 1, 0, 1, 0, 1, 0]
+loop_lengths = [1, 1]
+n_beam_splitters = len(loop_lengths) * len(input_state) - sum(loop_lengths)
+bs_angles = np.linspace(np.pi / 3, np.pi / 6, n_beam_splitters)
+ps_angles = np.linspace(np.pi / 3, np.pi / 5, n_beam_splitters)
+
+state = cudaq.get_state(TBI, bs_angles, ps_angles, input_state, loop_lengths)
+state.dump()
\ No newline at end of file
diff --git a/docs/sphinx/targets/python/photonics_tbi.py b/docs/sphinx/targets/python/photonics_tbi_sample.py
similarity index 71%
rename from docs/sphinx/targets/python/photonics_tbi.py
rename to docs/sphinx/targets/python/photonics_tbi_sample.py
index 712d34ead5..e381c6c46a 100644
--- a/docs/sphinx/targets/python/photonics_tbi.py
+++ b/docs/sphinx/targets/python/photonics_tbi_sample.py
@@ -14,24 +14,25 @@ def TBI(
     n_modes = len(input_state)
     level = sum(input_state) + 1  # qudit level
 
-    quds = [qudit(level) for _ in range(n_modes)]
+    qumodes = [qudit(level) for _ in range(n_modes)]
 
     for i in range(n_modes):
         for _ in range(input_state[i]):
-            plus(quds[i])
+            create(qumodes[i])
 
     counter = 0
     for j in loop_lengths:
         for i in range(n_modes - j):
-            beam_splitter(quds[i], quds[i + j], bs_angles[counter])
-            phase_shift(quds[i], ps_angles[counter])
+            beam_splitter(qumodes[i], qumodes[i + j], bs_angles[counter])
+            phase_shift(qumodes[i], ps_angles[counter])
             counter += 1
 
-    mz(quds)
+    mz(qumodes)
 
 
-input_state = [2, 1, 3, 1]
-loop_lengths = [1, 2]
+input_state = [1, 0, 1, 0, 1, 0, 1, 0]
+loop_lengths = [1, 1]
+
 n_beam_splitters = len(loop_lengths) * len(input_state) - sum(loop_lengths)
 bs_angles = np.linspace(np.pi / 3, np.pi / 6, n_beam_splitters)
 ps_angles = np.linspace(np.pi / 3, np.pi / 5, n_beam_splitters)
@@ -43,6 +44,3 @@ def TBI(
                       loop_lengths,
                       shots_count=1000000)
 counts.dump()
-
-state = cudaq.get_state(TBI, bs_angles, ps_angles, input_state, loop_lengths)
-state.dump()
diff --git a/python/cudaq/handlers/photonics_kernel.py b/python/cudaq/handlers/photonics_kernel.py
index 23d56bacfa..62a5b37caf 100644
--- a/python/cudaq/handlers/photonics_kernel.py
+++ b/python/cudaq/handlers/photonics_kernel.py
@@ -120,6 +120,40 @@ def qudit(level: int) -> PyQudit:
     return QuditManager.allocate(level)
 
 
+def create(qudit: PyQudit):
+    """
+    Apply create gate on the input qudit.
+    U|0> -> |1>, U|1> -> |2>, ..., and U|d> -> |d>
+
+    Args:
+        qudit: An instance of `PyQudit` class.
+
+    Raises:
+        RuntimeError: If the qudit level is not set.
+        Exception: If input argument is not instance of `PyQudit` class.
+    """
+    _check_args(qudit)
+    cudaq_runtime.photonics.apply_operation("create", [],
+                                            [[qudit.level, qudit.id]])
+
+
+def annihilate(qudit: PyQudit):
+    """
+    Apply annihilate gate on the input qudit.
+    U|0> -> |0>, U|1> -> |0>, ..., and U|d> -> |d-1>
+
+    Args:
+        qudit: An instance of `PyQudit` class.
+
+    Raises:
+        RuntimeError: If the qudit level is not set.
+        Exception: If input argument is not instance of `PyQudit` class.
+    """
+    _check_args(qudit)
+    cudaq_runtime.photonics.apply_operation("annihilate", [],
+                                            [[qudit.level, qudit.id]])
+
+
 def plus(qudit: PyQudit):
     """
     Apply plus gate on the input qudit. 
@@ -133,7 +167,7 @@ def plus(qudit: PyQudit):
         Exception: If input argument is not instance of `PyQudit` class.
     """
     _check_args(qudit)
-    cudaq_runtime.photonics.apply_operation('plusGate', [],
+    cudaq_runtime.photonics.apply_operation("plus", [],
                                             [[qudit.level, qudit.id]])
 
 
@@ -151,7 +185,7 @@ def phase_shift(qudit: PyQudit, phi: float):
         Exception: If input argument is not instance of `PyQudit` class.
     """
     _check_args(qudit)
-    cudaq_runtime.photonics.apply_operation('phaseShiftGate', [phi],
+    cudaq_runtime.photonics.apply_operation("phase_shift", [phi],
                                             [[qudit.level, qudit.id]])
 
 
@@ -170,7 +204,7 @@ def beam_splitter(q: PyQudit, r: PyQudit, theta: float):
         Exception: If input argument is not instance of `PyQudit` class.
     """
     _check_args([q, r])
-    cudaq_runtime.photonics.apply_operation('beamSplitterGate', [theta],
+    cudaq_runtime.photonics.apply_operation("beam_splitter", [theta],
                                             [[q.level, q.id], [r.level, r.id]])
 
 
@@ -209,7 +243,8 @@ class PhotonicsHandler(object):
     `qudit(level=N)` or a list of qudits. The qudits within a kernel must be of
     the same level.
 
-    Allowed quantum operations are: `plus`, `phase_shift`, `beam_splitter`, and `mz`.
+    Allowed quantum operations are: `create`, `annihilate`, `plus`,
+    `phase_shift`, `beam_splitter`, and `mz`.
     """
 
     def __init__(self, function):
@@ -221,11 +256,13 @@ def __init__(self, function):
         QuditManager.reset()
         self.kernelFunction = function
 
-        self.kernelFunction.__globals__['qudit'] = qudit
-        self.kernelFunction.__globals__['plus'] = plus
-        self.kernelFunction.__globals__['phase_shift'] = phase_shift
-        self.kernelFunction.__globals__['beam_splitter'] = beam_splitter
-        self.kernelFunction.__globals__['mz'] = mz
+        self.kernelFunction.__globals__["qudit"] = qudit
+        self.kernelFunction.__globals__["create"] = create
+        self.kernelFunction.__globals__["annihilate"] = annihilate
+        self.kernelFunction.__globals__["plus"] = plus
+        self.kernelFunction.__globals__["phase_shift"] = phase_shift
+        self.kernelFunction.__globals__["beam_splitter"] = beam_splitter
+        self.kernelFunction.__globals__["mz"] = mz
 
     def __call__(self, *args):
         with QuditManager():
diff --git a/python/tests/handlers/test_photonics_kernel.py b/python/tests/handlers/test_photonics_kernel.py
index 028b1af752..b61b258928 100644
--- a/python/tests/handlers/test_photonics_kernel.py
+++ b/python/tests/handlers/test_photonics_kernel.py
@@ -26,9 +26,9 @@ def test_qudit():
     @cudaq.kernel
     def kernel():
         q = qudit(level=4)
-        plus(q)
-        plus(q)
-        plus(q)
+        create(q)
+        create(q)
+        create(q)
         mz(q)
 
     counts = cudaq.sample(kernel)
@@ -44,10 +44,10 @@ def test_qudit_list():
 
     @cudaq.kernel
     def kernel():
-        qutrits = [qudit(3) for _ in range(2)]
-        plus(qutrits[0])
-        plus(qutrits[1])
-        mz(qutrits)
+        qumodes = [qudit(3) for _ in range(2)]
+        create(qumodes[0])
+        create(qumodes[1])
+        mz(qumodes)
 
     counts = cudaq.sample(kernel)
     assert len(counts) == 1
@@ -59,7 +59,7 @@ def test_qudit_invalid():
     @cudaq.kernel
     def kernel():
         q = [i for i in range(2)]
-        plus(q[0])
+        create(q[0])
 
     with pytest.raises(RuntimeError) as e:
         cudaq.sample(kernel)
@@ -70,16 +70,16 @@ def test_supported_gates():
 
     @cudaq.kernel
     def kernel():
-        quds = [qudit(5) for _ in range(3)]
+        qumodes = [qudit(5) for _ in range(3)]
 
-        plus(quds[0])
-        plus(quds[1])
-        plus(quds[2])
+        create(qumodes[0])
+        create(qumodes[1])
+        create(qumodes[2])
 
-        phase_shift(quds[1], 0.5)
-        beam_splitter(quds[0], quds[1], 1.3)
+        phase_shift(qumodes[1], 0.5)
+        beam_splitter(qumodes[0], qumodes[1], 1.3)
 
-        mz(quds)
+        mz(qumodes)
 
     counts = cudaq.sample(kernel)
     counts.dump()
@@ -92,7 +92,7 @@ def test_kernel_with_args():
     @cudaq.kernel
     def kernel_1f(theta: float):
         q = qudit(4)
-        plus(q)
+        create(q)
         phase_shift(q, theta)
         mz(q)
 
@@ -104,11 +104,11 @@ def kernel_1f(theta: float):
 
     @cudaq.kernel
     def kernel_2f(theta: float, phi: float):
-        quds = [qudit(3) for _ in range(2)]
-        plus(quds[0])
-        phase_shift(quds[0], theta)
-        beam_splitter(quds[0], quds[1], phi)
-        mz(quds)
+        qumodes = [qudit(3) for _ in range(2)]
+        create(qumodes[0])
+        phase_shift(qumodes[0], theta)
+        beam_splitter(qumodes[0], qumodes[1], phi)
+        mz(qumodes)
 
     result = cudaq.sample(kernel_2f, 0.7854, 0.3927)
     result.dump()
@@ -118,11 +118,11 @@ def kernel_2f(theta: float, phi: float):
 
     @cudaq.kernel
     def kernel_list(angles: List[float]):
-        quds = [qudit(2) for _ in range(3)]
-        plus(quds[0])
-        phase_shift(quds[1], angles[0])
-        phase_shift(quds[2], angles[1])
-        mz(quds)
+        qumodes = [qudit(2) for _ in range(3)]
+        create(qumodes[0])
+        phase_shift(qumodes[1], angles[0])
+        phase_shift(qumodes[2], angles[1])
+        mz(qumodes)
 
     result = cudaq.sample(kernel_list, [0.5236, 1.0472])
     result.dump()
@@ -136,7 +136,7 @@ def test_target_change():
     @cudaq.kernel
     def kernel():
         q = qudit(level=2)
-        plus(q)
+        create(q)
         mz(q)
 
     res = cudaq.sample(kernel)
@@ -185,7 +185,7 @@ def test_unsupported_types():
     @cudaq.kernel
     def kernel1():
         q = cudaq.qubit()
-        plus(q)
+        create(q)
 
     with pytest.raises(RuntimeError) as e:
         cudaq.sample(kernel1)
@@ -194,7 +194,7 @@ def kernel1():
     @cudaq.kernel
     def kernel2():
         q = cudaq.qvector(2)
-        plus(q[0])
+        create(q[0])
 
     with pytest.raises(Exception) as e:
         cudaq.sample(kernel2)
@@ -207,7 +207,7 @@ def test_target_handler():
     @cudaq.kernel
     def kernel():
         q = qudit(level=3)
-        plus(q)
+        create(q)
         mz(q)
 
     with pytest.raises(RuntimeError):
diff --git a/runtime/cudaq/qis/managers/photonics/PhotonicsExecutionManager.cpp b/runtime/cudaq/qis/managers/photonics/PhotonicsExecutionManager.cpp
index fc34484e5f..0e169461e7 100644
--- a/runtime/cudaq/qis/managers/photonics/PhotonicsExecutionManager.cpp
+++ b/runtime/cudaq/qis/managers/photonics/PhotonicsExecutionManager.cpp
@@ -131,7 +131,7 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
     state = qpp::kron(state, zeroState);
   }
 
-  /// @brief Allocate a set of `qudits` with a single call.
+  /// @brief Allocate a set of `qudits` (`qumodes`) with a single call.
   void allocateQudits(const std::vector<cudaq::QuditInfo> &qudits) override {
     for (auto &q : qudits)
       allocateQudit(q);
@@ -151,7 +151,7 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
   /// @brief Qudit deallocation method
   void deallocateQudit(const cudaq::QuditInfo &q) override {}
 
-  /// @brief Deallocate a set of `qudits` with a single call.
+  /// @brief Deallocate a set of `qudits` (`qumodes`) with a single call.
   void deallocateQudits(const std::vector<cudaq::QuditInfo> &qudits) override {}
 
   /// @brief Handler for when the photonics execution context changes
@@ -343,7 +343,33 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
 public:
   PhotonicsExecutionManager() {
 
-    instructions.emplace("plusGate", [&](const Instruction &inst) {
+    instructions.emplace("create", [&](const Instruction &inst) {
+      auto &[gateName, params, controls, qudits, spin_op] = inst;
+      auto target = qudits[0];
+      int d = target.levels;
+      qpp::cmat u{qpp::cmat::Zero(d, d)};
+      u(d - 1, d - 1) = 1;
+      for (int i = 1; i < d; i++) {
+        u(i, i - 1) = 1;
+      }
+      cudaq::info("Applying create on {}<{}>", target.id, target.levels);
+      state = qpp::apply(state, u, {target.id}, target.levels);
+    });
+
+    instructions.emplace("annihilate", [&](const Instruction &inst) {
+      auto &[gateName, params, controls, qudits, spin_op] = inst;
+      auto target = qudits[0];
+      int d = target.levels;
+      qpp::cmat u{qpp::cmat::Zero(d, d)};
+      u(0, 0) = 1;
+      for (int i = 0; i < d - 1; i++) {
+        u(i, i + 1) = 1;
+      }
+      cudaq::info("Applying annihilate on {}<{}>", target.id, target.levels);
+      state = qpp::apply(state, u, {target.id}, target.levels);
+    });
+
+    instructions.emplace("plus", [&](const Instruction &inst) {
       auto &[gateName, params, controls, qudits, spin_op] = inst;
       auto target = qudits[0];
       int d = target.levels;
@@ -352,11 +378,11 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
       for (int i = 1; i < d; i++) {
         u(i, i - 1) = 1;
       }
-      cudaq::info("Applying plusGate on {}<{}>", target.id, target.levels);
+      cudaq::info("Applying plus on {}<{}>", target.id, target.levels);
       state = qpp::apply(state, u, {target.id}, target.levels);
     });
 
-    instructions.emplace("beamSplitterGate", [&](const Instruction &inst) {
+    instructions.emplace("beam_splitter", [&](const Instruction &inst) {
       auto &[gateName, params, controls, qudits, spin_op] = inst;
       auto target1 = qudits[0];
       auto target2 = qudits[1];
@@ -364,12 +390,12 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
       const double theta = params[0];
       qpp::cmat BS{qpp::cmat::Zero(d * d, d * d)};
       beam_splitter(theta, BS);
-      cudaq::info("Applying beamSplitterGate on {}<{}> and {}<{}>", target1.id,
+      cudaq::info("Applying beam_splitter on {}<{}> and {}<{}>", target1.id,
                   target1.levels, target2.id, target2.levels);
       state = qpp::apply(state, BS, {target1.id, target2.id}, d);
     });
 
-    instructions.emplace("phaseShiftGate", [&](const Instruction &inst) {
+    instructions.emplace("phase_shift", [&](const Instruction &inst) {
       auto &[gateName, params, controls, qudits, spin_op] = inst;
       auto target = qudits[0];
       size_t d = target.levels;
@@ -379,8 +405,7 @@ class PhotonicsExecutionManager : public cudaq::BasicExecutionManager {
       for (size_t n = 0; n < d; n++) {
         PS(n, n) = std::exp(n * phi * i);
       }
-      cudaq::info("Applying phaseShiftGate on {}<{}>", target.id,
-                  target.levels);
+      cudaq::info("Applying phase_shift on {}<{}>", target.id, target.levels);
       state = qpp::apply(state, PS, {target.id}, target.levels);
     });
   }
diff --git a/runtime/cudaq/qis/managers/photonics/photonics_qis.h b/runtime/cudaq/qis/managers/photonics/photonics_qis.h
index ab4985e09b..ba2f04a095 100644
--- a/runtime/cudaq/qis/managers/photonics/photonics_qis.h
+++ b/runtime/cudaq/qis/managers/photonics/photonics_qis.h
@@ -15,18 +15,31 @@
 #include <vector>
 
 namespace cudaq {
+/// @brief The `create` gate
+// U|0> -> |1>, U|1> -> |2>, ..., and U|d> -> |d>
+template <std::size_t Levels>
+void create(qudit<Levels> &q) {
+  getExecutionManager()->apply("create", {}, {}, {{q.n_levels(), q.id()}});
+}
+
+/// @brief The `annihilate` gate
+// U|0> -> |0>, U|1> -> |0>, ..., and U|d> -> |d-1>
+template <std::size_t Levels>
+void annihilate(qudit<Levels> &q) {
+  getExecutionManager()->apply("annihilate", {}, {}, {{q.n_levels(), q.id()}});
+}
+
 /// @brief The `plus` gate
 // U|0> -> |1>, U|1> -> |2>, ..., and U|d> -> |0>
 template <std::size_t Levels>
 void plus(cudaq::qudit<Levels> &q) {
-  cudaq::getExecutionManager()->apply("plusGate", {}, {},
-                                      {{q.n_levels(), q.id()}});
+  cudaq::getExecutionManager()->apply("plus", {}, {}, {{q.n_levels(), q.id()}});
 }
 
 /// @brief The `phase shift` gate
 template <std::size_t Levels>
 void phase_shift(cudaq::qudit<Levels> &q, const double &phi) {
-  cudaq::getExecutionManager()->apply("phaseShiftGate", {phi}, {},
+  cudaq::getExecutionManager()->apply("phase_shift", {phi}, {},
                                       {{q.n_levels(), q.id()}});
 }
 
@@ -35,7 +48,7 @@ template <std::size_t Levels>
 void beam_splitter(cudaq::qudit<Levels> &q, cudaq::qudit<Levels> &r,
                    const double &theta) {
   cudaq::getExecutionManager()->apply(
-      "beamSplitterGate", {theta}, {},
+      "beam_splitter", {theta}, {},
       {{q.n_levels(), q.id()}, {r.n_levels(), r.id()}});
 }
 
diff --git a/runtime/cudaq/qis/state.cpp b/runtime/cudaq/qis/state.cpp
index 0253226265..461607f389 100644
--- a/runtime/cudaq/qis/state.cpp
+++ b/runtime/cudaq/qis/state.cpp
@@ -55,21 +55,35 @@ state::operator()(const std::initializer_list<std::size_t> &indices,
 
 std::complex<double> state::operator[](std::size_t idx) const {
   std::size_t numQubits = internal->getNumQubits();
+  std::size_t numElements = internal->getNumElements();
+
   if (!internal->isArrayLike()) {
     // Use amplitude accessor if linear indexing is not supported, e.g., tensor
     // network state.
     std::vector<int> basisState(numQubits, 0);
-    for (std::size_t i = 0; i < numQubits; ++i) {
-      if (idx & (1ULL << i))
-        basisState[(numQubits - 1) - i] = 1;
+    // Are we dealing with qudits or qubits?
+    if (std::log2(numElements) / numQubits > 1) {
+      for (std::size_t i = 0; i < numQubits; ++i) {
+        basisState[i] = 1; // TODO: This is a placeholder. We need to figure out
+                           // how to handle qudits.
+      }
+    } else {
+      for (std::size_t i = 0; i < numQubits; ++i) {
+        if (idx & (1ULL << i))
+          basisState[(numQubits - 1) - i] = 1;
+      }
     }
     return internal->getAmplitude(basisState);
   }
 
   std::size_t newIdx = 0;
-  for (std::size_t i = 0; i < numQubits; ++i)
-    if (idx & (1ULL << i))
-      newIdx |= (1ULL << ((numQubits - 1) - i));
+  if (std::log2(numElements) / numQubits > 1) {
+    newIdx = idx;
+  } else {
+    for (std::size_t i = 0; i < numQubits; ++i)
+      if (idx & (1ULL << i))
+        newIdx |= (1ULL << ((numQubits - 1) - i));
+  }
   return operator()({newIdx}, 0);
 }
 
diff --git a/unittests/photonics/PhotonicsTester.cpp b/unittests/photonics/PhotonicsTester.cpp
index 0fce4b70ad..552477f60e 100644
--- a/unittests/photonics/PhotonicsTester.cpp
+++ b/unittests/photonics/PhotonicsTester.cpp
@@ -15,21 +15,21 @@ TEST(PhotonicsTester, checkSimple) {
 
   struct test {
     auto operator()() __qpu__ {
-      cudaq::qvector<3> qutrits(2);
-      plus(qutrits[0]);
-      plus(qutrits[1]);
-      plus(qutrits[1]);
-      return mz(qutrits);
+      cudaq::qvector<3> qumodes(2);
+      create(qumodes[0]);
+      create(qumodes[1]);
+      create(qumodes[1]);
+      return mz(qumodes);
     }
   };
 
   struct test2 {
     void operator()() __qpu__ {
-      cudaq::qvector<3> qutrits(2);
-      plus(qutrits[0]);
-      plus(qutrits[1]);
-      plus(qutrits[1]);
-      mz(qutrits);
+      cudaq::qvector<3> qumodes(2);
+      create(qumodes[0]);
+      create(qumodes[1]);
+      create(qumodes[1]);
+      mz(qumodes);
     }
   };
 
@@ -51,15 +51,15 @@ TEST(PhotonicsTester, checkHOM) {
 
       constexpr std::array<std::size_t, 2> input_state{1, 1};
 
-      cudaq::qvector<3> quds(2); // |00>
+      cudaq::qvector<3> qumodes(2); // |00>
       for (std::size_t i = 0; i < 2; i++) {
         for (std::size_t j = 0; j < input_state[i]; j++) {
-          plus(quds[i]); // setting to  |11>
+          create(qumodes[i]); // setting to  |11>
         }
       }
 
-      beam_splitter(quds[0], quds[1], theta);
-      mz(quds);
+      beam_splitter(qumodes[0], qumodes[1], theta);
+      mz(qumodes);
     }
   };
 
@@ -93,18 +93,18 @@ TEST(PhotonicsTester, checkMZI) {
 
       constexpr std::array<std::size_t, 2> input_state{1, 0};
 
-      cudaq::qvector<3> quds(2); // |00>
+      cudaq::qvector<3> qumodes(2); // |00>
       for (std::size_t i = 0; i < 2; i++)
         for (std::size_t j = 0; j < input_state[i]; j++)
-          plus(quds[i]); // setting to  |10>
+          create(qumodes[i]); // setting to  |10>
 
-      beam_splitter(quds[0], quds[1], M_PI / 4);
-      phase_shift(quds[0], M_PI / 3);
+      beam_splitter(qumodes[0], qumodes[1], M_PI / 4);
+      phase_shift(qumodes[0], M_PI / 3);
 
-      beam_splitter(quds[0], quds[1], M_PI / 4);
-      phase_shift(quds[0], M_PI / 3);
+      beam_splitter(qumodes[0], qumodes[1], M_PI / 4);
+      phase_shift(qumodes[0], M_PI / 3);
 
-      mz(quds);
+      mz(qumodes);
     }
   };