diff --git a/tangelo/helpers/utils.py b/tangelo/helpers/utils.py
index 6ebc20324..ecb81f57b 100644
--- a/tangelo/helpers/utils.py
+++ b/tangelo/helpers/utils.py
@@ -73,9 +73,10 @@ def new_func(*args, **kwargs):
 
 
 # List all built-in backends supported
-all_backends = {"qulacs", "qiskit", "cirq", "braket", "projectq", "qdk", "pennylane"}
+all_backends = {"qulacs", "qiskit", "cirq", "braket", "projectq", "qdk", "pennylane", "sympy"}
 all_backends_simulator = {"qulacs", "qiskit", "cirq", "qdk"}
 sv_backends_simulator = {"qulacs", "qiskit", "cirq"}
+symbolic_backends = {"sympy"}
 
 # Dictionary mapping package names to their identifier in this module
 packages = {p: p for p in all_backends}
diff --git a/tangelo/linq/target/__init__.py b/tangelo/linq/target/__init__.py
index 91fc9c6be..64463fc32 100644
--- a/tangelo/linq/target/__init__.py
+++ b/tangelo/linq/target/__init__.py
@@ -17,10 +17,11 @@
 from .target_qiskit import QiskitSimulator
 from .target_qulacs import QulacsSimulator
 from .target_qdk import QDKSimulator
+from .target_sympy import SympySimulator
 from tangelo.helpers.utils import all_backends_simulator
 
 
-target_dict = {"qiskit": QiskitSimulator, "cirq": CirqSimulator, "qdk": QDKSimulator, "qulacs": QulacsSimulator}
+target_dict = {"qiskit": QiskitSimulator, "cirq": CirqSimulator, "qdk": QDKSimulator, "qulacs": QulacsSimulator, "sympy": SympySimulator}
 
 # Generate backend info dictionary
 backend_info = {sim_id: target_dict[sim_id].backend_info() for sim_id in all_backends_simulator}
diff --git a/tangelo/linq/target/target_sympy.py b/tangelo/linq/target/target_sympy.py
new file mode 100644
index 000000000..468ec569e
--- /dev/null
+++ b/tangelo/linq/target/target_sympy.py
@@ -0,0 +1,111 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# 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.
+
+import numpy as np
+
+from tangelo.linq import Circuit
+from tangelo.linq.target.backend import Backend
+from tangelo.linq.translator import translate_circuit as translate_c
+from tangelo.linq.translator import translate_operator
+
+
+class SympySimulator(Backend):
+
+    def __init__(self, n_shots=None, noise_model=None):
+        super().__init__(n_shots, noise_model)
+
+    def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None):
+        """This simulator manipulates symbolic expressions, i.e. gates can have
+        unspecified parameters (strings interpreted as variables). As with the
+        other simulators, it performs state preparation corresponding to the
+        input circuit, returns the frequencies of the different observables, and
+        either the statevector or None depending on if return_statevector is set
+        to True.
+
+        Args:
+            source_circuit (Circuit): A circuit in the abstract format to be
+                translated for the target backend.
+            return_statevector (bool): Option to return the statevector as well,
+                if available.
+            initial_statevector (array/matrix or sympy.physics.quantum.Qubit): A
+                valid statevector in the format supported by the target backend.
+
+        Returns:
+            dict: A dictionary mapping multi-qubit states to their corresponding
+                frequency.
+            sympy.Matrix: The symbolic statevector, if requested
+                by the user (if not, set to None).
+        """
+
+        from sympy import simplify
+        from sympy.physics.quantum import qapply
+        from sympy.physics.quantum.qubit import Qubit, matrix_to_qubit, \
+            qubit_to_matrix, measure_all
+
+        translated_circuit = translate_c(source_circuit, "sympy")
+
+        # Transform the initial_statevector if it is provided.
+        if initial_statevector is None:
+            python_statevector = Qubit("0"*(source_circuit.width))
+        elif isinstance(initial_statevector, Qubit):
+            python_statevector = initial_statevector
+        elif isinstance(initial_statevector, (np.ndarray, np.matrix)):
+            python_statevector = matrix_to_qubit(initial_statevector)
+        else:
+            raise ValueError(f"The {type(initial_statevector)} type for initial_statevector is not supported.")
+
+        # Deterministic circuit, run once.
+        state = qapply(translated_circuit * python_statevector)
+        self._current_state = state
+        python_statevector = qubit_to_matrix(state)
+
+        measurements = measure_all(state)
+
+        frequencies = dict()
+        for vec, prob in measurements:
+            prob = simplify(prob, tolerance=1e-4)
+            bistring = "".join(str(bit) for bit in reversed(vec.qubit_values))
+            frequencies[bistring] = prob
+
+        return (frequencies, python_statevector) if return_statevector else (frequencies, None)
+
+    def expectation_value_from_prepared_state(self, qubit_operator, n_qubits, prepared_state=None):
+        """Compute an expectation value using a representation of the state
+        using sympy functionalities.
+
+        Args:
+            qubit_operator (QubitOperator): a qubit operator in tangelo format
+            n_qubits (int): Number of qubits.
+            prepared_state (array/matrix or sympy.physics.quantum.Qubit): A
+                numpy or a sympy object representing the state. Internally, a
+                numpy object is transformed into the sympy representation.
+                Default is None, in this case it is set to the current state in
+                the simulator object.
+
+        Returns:
+            sympy.core.add.Add: Eigenvalue represented as a symbolic sum.
+        """
+
+        from sympy import simplify
+        from sympy.physics.quantum import qapply, Dagger
+
+        prepared_state = self._current_state if prepared_state is None else prepared_state
+        operator = translate_operator(qubit_operator, source="tangelo", target="sympy", n_qubits=n_qubits)
+        eigenvalue = qapply(Dagger(prepared_state) * operator * prepared_state)
+
+        return simplify(eigenvalue)
+
+    @staticmethod
+    def backend_info():
+        return {"statevector_available": True, "statevector_order": "lsq_first", "noisy_simulation": False}
diff --git a/tangelo/linq/tests/test_symbolic_simulator.py b/tangelo/linq/tests/test_symbolic_simulator.py
new file mode 100644
index 000000000..79bb6d2e7
--- /dev/null
+++ b/tangelo/linq/tests/test_symbolic_simulator.py
@@ -0,0 +1,82 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# 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.
+
+"""A test class to check that the simulator class functionalities are behaving
+as expected for the symbolic backend.
+"""
+
+import unittest
+from math import pi
+
+from tangelo.helpers.utils import assert_freq_dict_almost_equal
+from tangelo.linq import Gate, Circuit
+from tangelo.helpers.utils import installed_backends
+from tangelo.linq.target.target_sympy import SympySimulator
+
+
+class TestSymbolicSimulate(unittest.TestCase):
+
+    @unittest.skipIf("sympy" not in installed_backends, "Test Skipped: Sympy backend not available \n")
+    def test_simple_simulate(self):
+        """Test simulate of a simple rotation gate with a symbolic parameter."""
+
+        from sympy import symbols, cos, sin
+
+        simple_circuit = Circuit([Gate("RY", 0, parameter="alpha")])
+        backend = SympySimulator()
+        probs, _ = backend.simulate(simple_circuit, return_statevector=False)
+
+        alpha = symbols("alpha", real=True)
+
+        self.assertDictEqual(probs, {"0": (cos(alpha/2))**2, "1": (sin(alpha/2))**2})
+
+    @unittest.skipIf("sympy" not in installed_backends, "Test Skipped: Sympy backend not available \n")
+    def test_simulate_with_control(self):
+        """Test simulate of a control rotation gate with a symbolic parameter."""
+
+        from sympy import symbols, cos, sin
+
+        backend = SympySimulator()
+
+        no_action_circuit = Circuit([Gate("CRY", 1, 0, parameter="alpha")])
+        no_action_probs, _ = backend.simulate(no_action_circuit, return_statevector=False)
+
+        self.assertDictEqual(no_action_probs, {"00": 1.})
+
+        action_circuit = Circuit([Gate("X", 0), Gate("CRY", 1, 0, parameter="alpha")])
+        action_probs, _ = backend.simulate(action_circuit, return_statevector=False)
+        alpha = symbols("alpha", real=True)
+
+        self.assertDictEqual(action_probs, {"10": (cos(alpha/2))**2, "11": (sin(alpha/2))**2})
+
+    @unittest.skipIf("sympy" not in installed_backends, "Test Skipped: Sympy backend not available \n")
+    def test_evaluate_bell_state(self):
+        """Test the numerical evaluation to a known state (Bell state)."""
+
+        backend = SympySimulator()
+
+        variable_bell_circuit = Circuit([Gate("RY", 0, parameter="alpha"), Gate("CNOT", 1, 0)])
+        variable_bell_probs, _ = backend.simulate(variable_bell_circuit, return_statevector=False)
+
+        # Replace alpha by pi/2.
+        numerical_bell_probs = {
+            bitstring: prob.subs(list(prob.free_symbols)[0], pi/2) for
+            bitstring, prob in variable_bell_probs.items()
+        }
+
+        assert_freq_dict_almost_equal(numerical_bell_probs, {"00": 0.5, "11": 0.5}, atol=1e-3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tangelo/linq/tests/test_translator_circuit.py b/tangelo/linq/tests/test_translator_circuit.py
index 84acda14c..b89ec1255 100644
--- a/tangelo/linq/tests/test_translator_circuit.py
+++ b/tangelo/linq/tests/test_translator_circuit.py
@@ -501,6 +501,23 @@ def circuit(ops):
         # Compare statevectors
         np.testing.assert_array_almost_equal(v1, reference_big_lsq, decimal=6)
 
+    @unittest.skipIf("sympy" not in installed_backends, "Test Skipped: Sympy backend not available \n")
+    def test_to_sympy(self):
+        """Translate abtract format to sympy format."""
+
+        from sympy.physics.quantum.gate import HadamardGate, XGate, YGate, ZGate, CNotGate
+
+        # Equivalent native sympy circuit.
+        ref_circ = ZGate(3) * YGate(1) * XGate(0) * CNotGate(0, 1) * HadamardGate(2)
+
+        gates = [Gate("H", 2), Gate("CNOT", 1, control=0), Gate("X", 0), Gate("Y", 1), Gate("Z", 3)]
+        abs_circ = Circuit(gates)
+
+        # Generate the sympy circuit by translating from the abstract one.
+        translated_circuit = translate_c(abs_circ, "sympy")
+
+        self.assertEqual(translated_circuit, ref_circ)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/linq/tests/test_translator_perf.py b/tangelo/linq/tests/test_translator_perf.py
index 5b8de2dfb..061ae1d97 100644
--- a/tangelo/linq/tests/test_translator_perf.py
+++ b/tangelo/linq/tests/test_translator_perf.py
@@ -22,7 +22,7 @@
 
 from tangelo.linq import Gate, Circuit
 from tangelo.toolboxes.operators import QubitOperator
-from tangelo.helpers.utils import installed_backends
+from tangelo.helpers.utils import symbolic_backends
 from tangelo.linq import translate_operator, translate_circuit
 from tangelo.linq.translator.translate_qubitop import FROM_TANGELO as FROM_TANGELO_OP
 from tangelo.linq.translator.translate_qubitop import TO_TANGELO as TO_TANGELO_OP
@@ -47,12 +47,15 @@
 class PerfTranslatorTest(unittest.TestCase):
 
     def test_perf_operator(self):
-        """ Performance test with a reasonable large input for operator """
+        """ Performance test with a reasonable large input for operator.
+        Symbolic backends are not included in this test.
+        """
 
         print(f'\n[Performance Test :: linq operator format conversion]')
         print(f'\tInput size: n_qubits={n_qubits_op}, n_terms={n_terms}\n')
 
-        for f in FROM_TANGELO_OP:
+        perf_backends = FROM_TANGELO_OP.keys() - symbolic_backends
+        for f in perf_backends:
             try:
                 tstart = time.time()
                 target_op = translate_operator(tangelo_op, source="tangelo", target=f)
@@ -69,12 +72,15 @@ def test_perf_operator(self):
                     continue
 
     def test_perf_circuit(self):
-        """ Performance test with a reasonable large input for quantum circuit """
+        """ Performance test with a reasonable large input for quantum circuit.
+        Symbolic backends are not included in this test.
+        """
 
         print(f'\n[Performance Test :: linq circuit format conversion]')
         print(f'\tInput size: n_qubits={tangelo_c.width}, n_gates={tangelo_c.size}\n')
 
-        for f in FROM_TANGELO_C:
+        perf_backends = FROM_TANGELO_C.keys() - symbolic_backends
+        for f in perf_backends:
             try:
                 tstart = time.time()
                 target_c = translate_circuit(tangelo_c, source="tangelo", target=f)
diff --git a/tangelo/linq/translator/__init__.py b/tangelo/linq/translator/__init__.py
index a0fe09682..d3f1bc7cd 100644
--- a/tangelo/linq/translator/__init__.py
+++ b/tangelo/linq/translator/__init__.py
@@ -25,6 +25,7 @@
 from .translate_qubitop import translate_operator
 from .translate_circuit import translate_circuit
 from .translate_pennylane import get_pennylane_gates
+from .translate_sympy import get_sympy_gates
 
 
 def get_supported_gates():
@@ -40,5 +41,6 @@ def get_supported_gates():
     supported_gates["cirq"] = sorted(get_cirq_gates().keys())
     supported_gates["braket"] = sorted(get_braket_gates().keys())
     supported_gates["pennylane"] = sorted(get_pennylane_gates().keys())
+    supported_gates["sympy"] = sorted(get_sympy_gates().keys())
 
     return supported_gates
diff --git a/tangelo/linq/translator/translate_circuit.py b/tangelo/linq/translator/translate_circuit.py
index 65b9f4a54..74a0492c9 100644
--- a/tangelo/linq/translator/translate_circuit.py
+++ b/tangelo/linq/translator/translate_circuit.py
@@ -23,6 +23,7 @@
 from tangelo.linq.translator.translate_qiskit import translate_c_to_qiskit, translate_c_from_qiskit
 from tangelo.linq.translator.translate_qulacs import translate_c_to_qulacs
 from tangelo.linq.translator.translate_pennylane import translate_c_to_pennylane
+from tangelo.linq.translator.translate_sympy import translate_c_to_sympy
 
 
 FROM_TANGELO = {
@@ -34,7 +35,8 @@
     "qdk": translate_c_to_qsharp,
     "qiskit": translate_c_to_qiskit,
     "qulacs": translate_c_to_qulacs,
-    "pennylane": translate_c_to_pennylane
+    "pennylane": translate_c_to_pennylane,
+    "sympy": translate_c_to_sympy
 }
 
 TO_TANGELO = {
diff --git a/tangelo/linq/translator/translate_qubitop.py b/tangelo/linq/translator/translate_qubitop.py
index abe54d651..679ecd105 100644
--- a/tangelo/linq/translator/translate_qubitop.py
+++ b/tangelo/linq/translator/translate_qubitop.py
@@ -20,6 +20,7 @@
 from tangelo.linq.translator.translate_qulacs import translate_op_from_qulacs, translate_op_to_qulacs
 from tangelo.linq.translator.translate_pennylane import translate_op_from_pennylane, translate_op_to_pennylane
 from tangelo.linq.translator.translate_projectq import translate_op_from_projectq, translate_op_to_projectq
+from tangelo.linq.translator.translate_sympy import translate_op_to_sympy
 
 
 FROM_TANGELO = {
@@ -27,7 +28,8 @@
     "cirq": translate_op_to_cirq,
     "qulacs": translate_op_to_qulacs,
     "pennylane": translate_op_to_pennylane,
-    "projectq": translate_op_to_projectq
+    "projectq": translate_op_to_projectq,
+    "sympy": translate_op_to_sympy
 }
 
 TO_TANGELO = {
@@ -71,7 +73,7 @@ def translate_operator(qubit_operator, source, target, n_qubits=None):
             raise NotImplementedError(f"Qubit operator conversion from {source} to {target} is not supported.")
 
         # For translation functions that need an explicit number of qubits.
-        if target in {"qiskit"}:
+        if target in {"qiskit", "sympy"}:
             # The count_qubits function has no way to detect the number of
             # qubits when an operator is only a tensor product of I.
             if qubit_operator == QubitOperator((), qubit_operator.constant):
diff --git a/tangelo/linq/translator/translate_sympy.py b/tangelo/linq/translator/translate_sympy.py
new file mode 100644
index 000000000..79d2a1c25
--- /dev/null
+++ b/tangelo/linq/translator/translate_sympy.py
@@ -0,0 +1,230 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# 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.
+
+"""Functions helping with quantum circuit format conversion between abstract
+format and SYMPY format
+
+In order to produce an equivalent circuit for the target backend, it is necessary
+to account for:
+- how the gate names differ between the source backend to the target backend.
+- how the order and conventions for some of the inputs to the gate operations
+    may also differ.
+"""
+
+from tangelo.linq.helpers import pauli_of_to_string
+
+
+def rx_gate(target, theta):
+    """Implementation of the RX gate with a unitary matrix.
+
+    Args:
+        target (int): Target qubit.
+        theta (float): Rotation angle.
+
+    Returns:
+        sympy.physics.quantum.gate.UGate: Self-explanatory.
+    """
+
+    from sympy import ImmutableMatrix, sin, cos, I
+    from sympy.physics.quantum.gate import UGate
+
+    cos_term = cos(theta / 2)
+    sin_term = -I*sin(theta / 2)
+    rx_matrix = ImmutableMatrix([[cos_term, sin_term], [sin_term, cos_term]])
+
+    return UGate(target, rx_matrix)
+
+
+def ry_gate(target, theta):
+    """Implementation of the RY gate with a unitary matrix.
+
+    Args:
+        target (int): Target qubit.
+        theta (float): Rotation angle.
+
+    Returns:
+        sympy.physics.quantum.gate.UGate: Self-explanatory.
+    """
+
+    from sympy import ImmutableMatrix, sin, cos
+    from sympy.physics.quantum.gate import UGate
+
+    cos_term = cos(theta / 2)
+    sin_term = sin(theta / 2)
+    ry_matrix = ImmutableMatrix([[cos_term, -sin_term], [sin_term, cos_term]])
+
+    return UGate(target, ry_matrix)
+
+
+def rz_gate(target, theta):
+    """Implementation of the RZ gate with a unitary matrix.
+
+    Args:
+        target (int): Target qubit.
+        theta (float): Rotation angle.
+
+    Returns:
+        sympy.physics.quantum.gate.UGate: Self-explanatory.
+    """
+
+    from sympy import ImmutableMatrix, I, exp
+    from sympy.physics.quantum.gate import UGate
+
+    rz_matrix = ImmutableMatrix([[exp(- I * theta / 2), 0], [0, exp(I * theta / 2)]])
+
+    return UGate(target, rz_matrix)
+
+
+def p_gate(target, theta):
+    """Implementation of the PHASE gate with a unitary matrix.
+
+    Args:
+        target (int): Target qubit.
+        theta (float): Phase parameter.
+
+    Returns:
+        sympy.physics.quantum.gate.UGate: Self-explanatory.
+    """
+
+    from sympy import ImmutableMatrix, I, exp
+    from sympy.physics.quantum.gate import UGate
+
+    p_matrix = ImmutableMatrix([[1, 0], [0, exp(I * theta)]])
+
+    return UGate(target, p_matrix)
+
+
+def controlled_gate(gate_function):
+    """Change a one-qubit gate to a controlled gate.
+
+    Args:
+        gate_function (sympy.physics.quantum.gate): Class for a specific gate.
+
+    Returns:
+        function: Function wrapping a gate with CGate.
+    """
+
+    from sympy.physics.quantum.gate import CGate
+
+    def cgate(control, target, *args, **kwargs):
+        return CGate(control, gate_function(target, *args, **kwargs))
+
+    return cgate
+
+
+def get_sympy_gates():
+    """Map gate name of the abstract format to the equivalent methods of the
+    sympy backend (symbolic simulation).
+    """
+
+    import sympy.physics.quantum.gate as SYMPYGate
+
+    GATE_SYMPY = dict()
+    GATE_SYMPY["H"] = SYMPYGate.HadamardGate
+    GATE_SYMPY["X"] = SYMPYGate.XGate
+    GATE_SYMPY["Y"] = SYMPYGate.YGate
+    GATE_SYMPY["Z"] = SYMPYGate.ZGate
+    GATE_SYMPY["S"] = SYMPYGate.PhaseGate
+    GATE_SYMPY["T"] = SYMPYGate.TGate
+    GATE_SYMPY["PHASE"] = p_gate
+    GATE_SYMPY["SWAP"] = SYMPYGate.SwapGate
+    GATE_SYMPY["RX"] = rx_gate
+    GATE_SYMPY["RY"] = ry_gate
+    GATE_SYMPY["RZ"] = rz_gate
+    GATE_SYMPY["RX"] = rx_gate
+    GATE_SYMPY["CH"] = controlled_gate(SYMPYGate.HadamardGate)
+    GATE_SYMPY["CNOT"] = SYMPYGate.CNotGate
+    GATE_SYMPY["CX"] = SYMPYGate.CNotGate
+    GATE_SYMPY["CY"] = controlled_gate(SYMPYGate.YGate)
+    GATE_SYMPY["CZ"] = controlled_gate(SYMPYGate.ZGate)
+    GATE_SYMPY["CRX"] = controlled_gate(rx_gate)
+    GATE_SYMPY["CRY"] = controlled_gate(ry_gate)
+    GATE_SYMPY["CRZ"] = controlled_gate(rz_gate)
+    GATE_SYMPY["CS"] = controlled_gate(SYMPYGate.PhaseGate)
+    GATE_SYMPY["CT"] = controlled_gate(SYMPYGate.TGate)
+    GATE_SYMPY["CPHASE"] = controlled_gate(p_gate)
+
+    return GATE_SYMPY
+
+
+def translate_c_to_sympy(source_circuit):
+    """Take in an abstract circuit, return a quantum circuit object as defined
+    in the Python SYMPY SDK.
+
+    Args:
+        source_circuit: quantum circuit in the abstract format.
+
+    Returns:
+        SYMPY.circuits.Circuit: quantum circuit in Python SYMPY SDK format.
+    """
+
+    from sympy import symbols
+
+    GATE_SYMPY = get_sympy_gates()
+
+    # Identity as an empty circuit.
+    target_circuit = 1
+
+    # Map the gate information properly.
+    for gate in reversed(source_circuit._gates):
+        # If the parameter is a string, we use it as a variable.
+        if gate.parameter and isinstance(gate.parameter, str):
+            gate.parameter = symbols(gate.parameter, real=True)
+
+        if gate.name in {"H", "X", "Y", "Z"}:
+            target_circuit *= GATE_SYMPY[gate.name](gate.target[0])
+        elif gate.name in {"T", "S"} and gate.parameter == "":
+            target_circuit *= GATE_SYMPY[gate.name](gate.target[0])
+        elif gate.name in {"PHASE", "RX", "RY", "RZ"}:
+            target_circuit *= GATE_SYMPY[gate.name](gate.target[0], gate.parameter)
+        elif gate.name in {"CNOT", "CH", "CX", "CY", "CZ", "CS", "CT"}:
+            target_circuit *= GATE_SYMPY[gate.name](gate.control[0], gate.target[0])
+        elif gate.name in {"SWAP"}:
+            target_circuit *= GATE_SYMPY[gate.name](gate.target[0], gate.target[1])
+        elif gate.name in {"CRX", "CRY", "CRZ", "CPHASE"}:
+            target_circuit *= GATE_SYMPY[gate.name](gate.control[0], gate.target[0], gate.parameter)
+        else:
+            raise ValueError(f"Gate '{gate.name}' not supported on backend SYMPY")
+
+    return target_circuit
+
+
+def translate_op_to_sympy(qubit_operator, n_qubits):
+    """Helper function to translate a Tangelo QubitOperator to a sympy linear
+    combination of tensor products.
+
+    Args:
+        qubit_operator (tangelo.toolboxes.operators.QubitOperator): Self-explanatory.
+        n_qubits (int): The number of qubit the operator acts on.
+
+    Returns:
+        sympy.core.add.Add: Summation of sympy.physics.quantum.TensorProduct
+            objects.
+    """
+    from sympy import Identity
+    from sympy.physics.paulialgebra import Pauli
+    from sympy.physics.quantum import TensorProduct
+
+    # Pauli string to sympy Pauli algebra objects.
+    map_to_paulis = {"I": Identity(1), "X": Pauli(1), "Y": Pauli(2), "Z": Pauli(3)}
+
+    # Contruct the TensorProduct objects.
+    sum_tensor_paulis = 0.
+    for term_tuple, coeff in qubit_operator.terms.items():
+        term_string = pauli_of_to_string(term_tuple, n_qubits)
+        paulis = [map_to_paulis[p] for p in term_string[::-1]]
+        tensor_paulis = TensorProduct(*paulis)
+        sum_tensor_paulis += coeff * tensor_paulis
+
+    return sum_tensor_paulis