diff --git a/backend/ibex/core/data_manipulation_methods.py b/backend/ibex/core/data_manipulation_methods.py
index 70d89a9c..978b394d 100644
--- a/backend/ibex/core/data_manipulation_methods.py
+++ b/backend/ibex/core/data_manipulation_methods.py
@@ -221,3 +221,44 @@ class DataManipulationMethodsResponse(BaseModel):
 
 simple_data_operations_description.method_parameters.append(data_operations_parameter)
 available_methods.data_manipulation_methods.append(simple_data_operations_description)
+
+# ====================== SIGNAL DATA OPERATIONS ======================
+
+signal_data_operations_description = DataManipulationOperation(
+    name="Signal Data Operations",
+    description="Ordered list of signal operations applied to the dataset. "
+    "Execution order is determined by the order of parameters in the request. "
+    "Each value is a URI pointing to a signal.",
+    method_parameters=[],
+)
+
+signal_data_operations_parameter = DataManipulationParameter(
+    human_readable_name="Operations",
+    name="operations",
+    description="Ordered list of signal operations applied to every data point.",
+    type="list[object]",
+    group_label="Operation",
+    fields=[
+        DataManipulationParameter(
+            human_readable_name="Type",
+            name="operation_type",
+            description="Type of operation",
+            type="string",
+            possible_values=[
+                PossibleValue(value="add", description="Addition"),
+                PossibleValue(value="sub", description="Subtraction"),
+                PossibleValue(value="mul", description="Multiplication"),
+                PossibleValue(value="div", description="Division"),
+            ],
+        ),
+        DataManipulationParameter(
+            human_readable_name="Signal URI",
+            name="operation_value",
+            description="URI to signal data for the operation",
+            type="string",
+        ),
+    ],
+)
+
+signal_data_operations_description.method_parameters.append(signal_data_operations_parameter)
+available_methods.data_manipulation_methods.append(signal_data_operations_description)
diff --git a/backend/ibex/data_source/imas_python_source.py b/backend/ibex/data_source/imas_python_source.py
index 3017a913..928b9e8e 100644
--- a/backend/ibex/data_source/imas_python_source.py
+++ b/backend/ibex/data_source/imas_python_source.py
@@ -5,7 +5,6 @@
 import imas  # type: ignore
 import numpy as np  # type: ignore
 import re  # type: ignore
-from copy import copy  # type: ignore
 from idstools.database import DBMaster  # type: ignore
 from imas.ids_metadata import IDSMetadata  # type: ignore
 from imas.ids_primitive import (
@@ -52,6 +51,8 @@
     apply_gaussian_filter,
     apply_simple_operations,
     resolve_irregular_coordinate_data_shape,
+    apply_signal_operations,
+    combine_signal_units,
 )
 from ibex.core.data_manipulation_methods import SmoothingMethod, InterpolationMethod
 from ibex.endpoints.schemas.request_data_schemas import PlotDataRequestModel
@@ -1055,6 +1056,7 @@ def get_plot_data(self, plot_data_query: PlotDataRequestModel) -> dict:
                         }
                         coordinates_to_be_returned.append(c)
             first_value = find_first_value_in_list(ids_data)
+            result_unit = first_value.metadata.units or ""
             data_to_be_returned = convert_ids_data_into_numpy_array(ids_data)
 
             if first_value.metadata.ndim == 2:
@@ -1106,8 +1108,12 @@ def get_plot_data(self, plot_data_query: PlotDataRequestModel) -> dict:
 
             # ============= END data smoothing =============
 
-            # ============= BEGIN resample data onto new time vector =============
+            try:
+                original_data_shape = np.asarray(data_to_be_returned).shape
+            except ValueError:
+                original_data_shape = "irregular"
 
+            # ============= BEGIN resample data onto new time vector =============
             def convert_to_lists(data):
                 if isinstance(data, list):
                     return [convert_to_lists(d) for d in data]
@@ -1116,6 +1122,10 @@ def convert_to_lists(data):
                 else:
                     return data
 
+            # list of dicts used when combining signals after
+            # {uri:str, data:list[*], interpolated_data: list[*]}
+            others_signals_data = {}
+
             if plot_data_query.interpolate_over:
                 # check for non-interpolatable coordinates
                 if plot_data_query.interpolation_method != InterpolationMethod.EXACT_VALUE:
@@ -1153,6 +1163,10 @@ def convert_to_lists(data):
 
                 original_coord_values.reverse()
 
+                store_other_signals_data = False  # used for signal combining after interpolation
+                if plot_data_query.signal_operations:
+                    store_other_signals_data = True
+
                 for _uri in plot_data_query.interpolate_over:
                     _uri_obj = IMAS_URI(_uri)
 
@@ -1165,11 +1179,20 @@ def convert_to_lists(data):
                                 "IDS name and node path should be the same for source and target URI when interpolating data"
                             )
 
-                    new_plot_data_query = copy(plot_data_query)
-                    new_plot_data_query.uri = _uri
-                    new_plot_data_query.interpolate_over = None
-                    new_plot_data_query.smoothing_method = None
-                    interpolate_to_coordinates = self.get_plot_data(new_plot_data_query)["data"]["coordinates"]
+                    new_plot_data_query = PlotDataRequestModel(uri=_uri)
+                    # interpolate_to will be used later with signal combining
+                    interpolate_to = self.get_plot_data(new_plot_data_query)["data"]
+                    interpolate_to_coordinates = interpolate_to["coordinates"]
+                    if store_other_signals_data:
+                        others_signals_data[_uri] = {
+                            "uri": _uri,
+                            "data": pad_to_rectangular(interpolate_to["value"]),
+                            "coordinates": [
+                                sorted(set(flatten(convert_to_lists(c["value"])))) for c in interpolate_to_coordinates
+                            ],
+                            "shape": interpolate_to["shape"],
+                            "unit": interpolate_to["unit"],
+                        }
 
                     if len(interpolate_to_coordinates) != len(coordinates_to_be_returned):
                         message = "Interpolation error. Source and target nodes have different number of coordinates."
@@ -1215,10 +1238,112 @@ def convert_to_lists(data):
 
             # ============= END resample data onto new time vector =============
 
+            # ============= BEGIN signal operations =============
+            #
+            # Steps performed in this block:
+            # 1. Collect all signal URIs referenced in signal_operations
+            # 2. Pad data to rectangular if shape is irregular
+            # 3. For each signal URI, fetch and prepare data:
+            #    a. If the signal was already interpolated (stored during
+            #       interpolation phase), skip fetching
+            #    b. Otherwise fetch the signal and check shape compatibility
+            # 4. Prepare operand data for each signal:
+            #    a. If interpolation was requested, resample onto common coords
+            #    b. Otherwise normalize raw signal data to a NumPy array
+            # 5. Build a flat uri->array dict and apply all signal operations
+
+            if plot_data_query.signal_operations:
+                # Step 1: extract unique signal URIs from operation strings
+                signal_op_uris = set()
+                for op_str in plot_data_query.signal_operations or []:
+                    _, op_uri = op_str.split(":", 1)
+                    signal_op_uris.add(op_uri)
+
+                # Step 2: ensure rectangular data for downstream processing
+                if original_data_shape == "irregular":
+                    data_to_be_returned = pad_to_rectangular(data_to_be_returned)
+
+                # Step 3: fetch and prepare each signal referenced in operations
+                for signal_uri in signal_op_uris:
+                    if signal_uri not in others_signals_data:
+                        # Signal was not pre-loaded during interpolation phase.
+                        # Fetch it now and verify shape compatibility.
+                        request = PlotDataRequestModel(uri=signal_uri)
+                        other_signal = self.get_plot_data(request)
+
+                        # Comparing signal shapes
+                        if (
+                            other_signal["data"]["shape"] == "irregular"
+                            or other_signal["data"]["shape"] != original_data_shape
+                        ):
+                            msg = (
+                                f"Cannot apply operation on signal {signal_uri} without interpolation. "
+                                "Signal and data shapes do not match. "
+                                "Try interpolating the signal onto the data shape."
+                            )
+                            raise InvalidParametersException(msg)
+
+                        # Comparing coordinates
+                        coordinates_match = self._coordinates_match(
+                            coordinates_to_be_returned, other_signal["data"]["coordinates"]
+                        )
+
+                        if not coordinates_match:
+                            msg = (
+                                f"Cannot apply operation on signal {signal_uri} without interpolation. "
+                                "Coordinates do not match. "
+                                "Try interpolating the signal onto the data coordinates."
+                            )
+                            raise InvalidParametersException(msg)
+
+                        others_signals_data[signal_uri] = {
+                            "uri": request.uri,
+                            "data": other_signal["data"]["value"],
+                            "coordinates": [
+                                sorted(set(flatten(convert_to_lists(c["value"]))))
+                                for c in other_signal["data"]["coordinates"]
+                            ],
+                            "shape": other_signal["data"]["shape"],
+                            "unit": other_signal["data"]["unit"],
+                        }
+
+                    # Step 4: prepare operand data (resampled or raw)
+                    if plot_data_query.interpolate_over:
+                        # Resample signal data onto the common coordinate grid
+                        signal_data = resample_data_without_interpolation(
+                            tuple(reversed(others_signals_data[signal_uri]["coordinates"])),
+                            others_signals_data[signal_uri]["data"],
+                            tuple(common_coords_values),
+                        )
+                    else:
+                        signal_data = others_signals_data[signal_uri]["data"]
+
+                    others_signals_data[signal_uri]["data"] = np.asarray(signal_data)
+
+                # Step 5: flatten dict and apply operations in order
+                signal_data_by_uri = {uri: info["data"] for uri, info in others_signals_data.items()}
+                for operation in plot_data_query.signal_operations:
+                    operation_type, signal_uri = operation.split(":", 1)
+                    result_unit = combine_signal_units(
+                        result_unit,
+                        others_signals_data[signal_uri]["unit"],
+                        operation_type,
+                    )
+                data_to_be_returned = apply_signal_operations(
+                    data_to_be_returned, plot_data_query.signal_operations, signal_data_by_uri
+                )
+
+            # ============= END signal operations =============
+
+            # Interpolation and signal operations can change the data shape.
+            # The response's ``shape`` describes the processed data before
+            # downsampling, rather than the raw shape used for compatibility
+            # checks above.
             try:
-                original_data_shape = np.asarray(data_to_be_returned).shape
+                processed_data_shape = np.asarray(data_to_be_returned).shape
             except ValueError:
-                original_data_shape = "irregular"
+                processed_data_shape = "irregular"
+
             # Downsample only 1D data
             if first_value.metadata.ndim == 1:
                 if coordinates_to_be_returned[0]["target"].split("/")[-1] == f"{node_path.split('/')[-1]}":
@@ -1250,8 +1375,8 @@ def convert_to_lists(data):
             result = {
                 "data": {
                     "name": node_path.split("/")[-1],
-                    "unit": first_value.metadata.units,
-                    "shape": original_data_shape,
+                    "unit": result_unit,
+                    "shape": processed_data_shape,
                     "downsampled_shape": downsampled_shape,
                     "ndim": first_value.metadata.ndim,
                     "path": f"#{ids}/{node_path}",
@@ -1272,6 +1397,17 @@ def convert_to_lists(data):
                 coordinate["coordinates"] = new_shape_factors_list
         return result
 
+    def _coordinates_match(self, coordinates_1, coordinates_2):
+        if len(coordinates_1) != len(coordinates_2):
+            return False
+
+        coordinates_match = all(
+            coordinate_1["name"] == coordinate_2["name"]
+            and np.array_equal(np.asarray(coordinate_1["value"]), np.asarray(coordinate_2["value"]), equal_nan=True)
+            for coordinate_1, coordinate_2 in zip(coordinates_1, coordinates_2)
+        )
+        return coordinates_match
+
     def _is_empty(self, seq):
         """Checks if list is essentially empty (contains only empty lists or empty strings)"""
         if isinstance(seq, (IDSNumericArray, IDSString0D, IDSString1D, IDSComplex0D, IDSFloat0D, IDSInt0D)):
diff --git a/backend/ibex/data_source/imas_python_source_utils.py b/backend/ibex/data_source/imas_python_source_utils.py
index 5a48cc65..3572f715 100644
--- a/backend/ibex/data_source/imas_python_source_utils.py
+++ b/backend/ibex/data_source/imas_python_source_utils.py
@@ -183,18 +183,21 @@ def apply_gaussian_filter(data: list | np.ndarray, sigma, axis: int | None = Non
         raise InvalidParametersException(msg)
 
 
-def _safe_division(data, divisor):
-    if divisor == 0:
-        raise InvalidParametersException("Division by zero is not allowed")
-    return data / divisor
-
-
 def _safe_root(data, exponent):
     if exponent == 0:
         raise InvalidParametersException("Root by zero is not allowed")
     return np.power(data, 1 / exponent)
 
 
+def _safe_division(data, divisor):
+    if isinstance(divisor, np.ndarray):
+        if np.any(divisor == 0):
+            raise InvalidParametersException("Division by zero is not allowed")
+    elif divisor == 0:
+        raise InvalidParametersException("Division by zero is not allowed")
+    return data / divisor
+
+
 _SIMPLE_OPERATIONS_FUNCTIONS = {
     "add": _op.add,
     "sub": _op.sub,
@@ -204,6 +207,40 @@ def _safe_root(data, exponent):
     "root": _safe_root,
 }
 
+_SIGNAL_OPERATIONS_FUNCTIONS = {
+    "add": _op.add,
+    "sub": _op.sub,
+    "mul": _op.mul,
+    "div": _safe_division,
+}
+
+
+def combine_signal_units(left_unit: str, right_unit: str, operation: str) -> str:
+    """Return the unit produced by a binary signal operation."""
+    left_unit = left_unit or ""
+    right_unit = right_unit or ""
+
+    if operation in {"add", "sub"}:
+        if left_unit != right_unit:
+            raise InvalidParametersException(
+                f"Cannot {operation} signals with different units ({left_unit!r} and {right_unit!r})"
+            )
+        return left_unit
+    if operation == "mul":
+        if not left_unit:
+            return right_unit
+        if not right_unit:
+            return left_unit
+        return f"{left_unit}*{right_unit}"
+    if operation == "div":
+        if left_unit == right_unit:
+            return ""
+        if not right_unit:
+            return left_unit
+        return f"{left_unit}/{right_unit}"
+
+    raise InvalidParametersException(f"Unknown operation type: {operation}")
+
 
 def apply_simple_operations(data: list | np.ndarray, operations: list[str]):
     """
@@ -230,6 +267,34 @@ def apply_simple_operations(data: list | np.ndarray, operations: list[str]):
         raise InvalidParametersException(msg)
 
 
+def apply_signal_operations(data: list | np.ndarray, operations: list[str], signal_data_by_uri: dict):
+    """
+    Apply signal operations to data in the order given.
+    Each operation is a string in the format 'type:uri', e.g. 'add:some/imas/uri', 'mul:other/uri'.
+    :param data: Input data
+    :param operations: List of operations and URIs divided by colon (:)
+    :param signal_data_by_uri: Dict mapping signal URIs to their interpolated data arrays.
+    :return: Data after operation
+    """
+
+    if isinstance(data, list):
+        data = np.array(data)
+    if isinstance(data, (np.ndarray, IDSNumericArray)):
+        result = data
+        for op_str in operations:
+            op_type, uri = op_str.split(":", 1)
+            value = signal_data_by_uri[uri]
+            func = _SIGNAL_OPERATIONS_FUNCTIONS.get(op_type)
+            if func is None:
+                raise InvalidParametersException(f"Unknown operation type: {op_type}")
+
+            result = func(result, value)
+        return result
+    else:
+        msg = "Signal operations can be executed only on numeric arrays, not single values or strings."
+        raise InvalidParametersException(msg)
+
+
 def union_arrays(data: list):
     return reduce(np.union1d, data)
 
@@ -253,8 +318,8 @@ def calculate_coordinate_shapes(shape: list[int], dims: int):
     :return: List of shapes for each coordinate.
     """
 
-    if dims < 0 or dims >= len(shape):
-        raise ValueError("dims must be >= 0 and < len(shape)")
+    if dims < 0 or dims > len(shape):
+        raise ValueError("dims must be >= 0 and <= len(shape)")
 
     # Base dimensions (dimensions added by AoS in path), e.g. [4, 5]
     base = shape[:dims]
@@ -399,16 +464,7 @@ def resample_data_without_interpolation(original_coords, data, target_coords):
     target_indices = []
 
     for orig, target in zip(original_coords, target_coords):
-        # Build dictionary: coordinate -> target index
-        lookup = {}
-        for i, value in enumerate(target):
-            lookup[value] = i
-
-        axis_indices = []
-
-        for value in orig:
-            axis_indices.append(lookup[value])
-
+        axis_indices = np.searchsorted(np.asarray(target), np.asarray(orig))
         target_indices.append(axis_indices)
 
     # Create mesh
diff --git a/backend/ibex/endpoints/schemas/request_data_schemas.py b/backend/ibex/endpoints/schemas/request_data_schemas.py
index 3acbc2f7..f6eb0475 100644
--- a/backend/ibex/endpoints/schemas/request_data_schemas.py
+++ b/backend/ibex/endpoints/schemas/request_data_schemas.py
@@ -1,7 +1,7 @@
 from pydantic import BaseModel, Field, model_validator
 from typing import Optional, List
 from ibex.core.data_manipulation_methods import available_methods, SmoothingMethod
-from ibex.data_source.imas_python_source_utils import _SIMPLE_OPERATIONS_FUNCTIONS
+from ibex.data_source.imas_python_source_utils import _SIMPLE_OPERATIONS_FUNCTIONS, _SIGNAL_OPERATIONS_FUNCTIONS
 from enum import Enum
 
 
@@ -97,6 +97,10 @@ class PlotDataBasicParameters(BaseModel):
         default=None,
         description="Ordered list of scalar operations in format 'type:value' e.g. 'add:10'",
     )
+    signal_operations: Optional[List[str]] = Field(
+        default=None,
+        description="Ordered list of signal operations in format 'type:uri' e.g. 'add:some/imas/uri'",
+    )
 
 
 class PlotDataRequestModel(
@@ -122,6 +126,22 @@ def validate_operations_format(self) -> "PlotDataRequestModel":
                     raise ValueError(f"Invalid operation value: '{value_str}' in '{op}'. Value must be a number.")
         return self
 
+    @model_validator(mode="after")
+    def validate_signal_operations_format(self) -> "PlotDataRequestModel":
+        if self.signal_operations:
+            valid_operations = set(_SIGNAL_OPERATIONS_FUNCTIONS.keys())
+            for op in self.signal_operations:
+                if ":" not in op:
+                    raise ValueError(f"Invalid signal operation format: '{op}'. Expected 'operation:uri'")
+                op_type, uri = op.split(":", 1)
+                if op_type not in valid_operations:
+                    raise ValueError(
+                        f"Unknown signal operation type: '{op_type}'. Valid types: {', '.join(sorted(valid_operations))}"
+                    )
+                if not uri.strip():
+                    raise ValueError(f"Invalid signal operation: '{op}'. URI must not be empty.")
+        return self
+
     @model_validator(mode="after")
     def validate_gaussian_smoothing_parameters(self) -> "PlotDataRequestModel":
         if self.smoothing_method == SmoothingMethod.GAUSSIAN_FILTER and self.gaussian_smoothing_sigma is None:
diff --git a/backend/pyproject.toml b/backend/pyproject.toml
index bb02a7b4..715d9f85 100644
--- a/backend/pyproject.toml
+++ b/backend/pyproject.toml
@@ -57,9 +57,9 @@ dependencies = [
 [project.optional-dependencies]
 # these self-dependencies are available since pip 21.2
 all = [
-    "ibex[test]",
-    "ibex[linting]",
-    "ibex[docs]",
+    "imas-ibex[test]",
+    "imas-ibex[linting]",
+    "imas-ibex[docs]",
 ]
 
 test = [
diff --git a/backend/tests/conftest.py b/backend/tests/conftest.py
index 6f630722..1bfebd51 100644
--- a/backend/tests/conftest.py
+++ b/backend/tests/conftest.py
@@ -11,12 +11,15 @@
 @pytest.fixture(scope="session")
 def interpolation_entry_path_directory(tmp_path_factory):
     tmp_path = tmp_path_factory.mktemp("interpolation_testdb")
+    rand_generator = np.random.default_rng(2137)
 
     with imas.DBEntry(f"imas:hdf5?path={tmp_path}/interpolation_db_1", mode="w") as entry:
         eq = entry.factory.equilibrium()
 
         eq.ids_properties.homogeneous_time = 1
         eq.time = np.asarray([1, 2, 3, 4], dtype=float)
+        eq.vacuum_toroidal_field.r0 = 1.0
+        eq.vacuum_toroidal_field.b0 = np.asarray([0.1, 0.2, 0.3, 0.4], dtype=float)
         eq.time_slice.resize(4)
         for ts in eq.time_slice:
             ts.profiles_1d.psi = np.asarray([1.0, 1.2, 1.3, 1.4, 1.5, 1.6])
@@ -29,7 +32,7 @@ def interpolation_entry_path_directory(tmp_path_factory):
 
             ts.profiles_2d.resize(2)
             for p2d in ts.profiles_2d:
-                p2d.psi = np.asarray(np.random.rand(3, 3))
+                p2d.psi = np.asarray(rand_generator.random((3, 3)))
                 p2d.grid.dim1 = np.asarray([1.0, 2.0, 3.0], dtype=float)
                 p2d.grid.dim2 = np.asarray([1.0, 2.0, 3.0], dtype=float)
         entry.put(eq)
@@ -65,6 +68,8 @@ def interpolation_entry_path_directory(tmp_path_factory):
 
         eq.ids_properties.homogeneous_time = 1
         eq.time = np.asarray([1, 2, 3], dtype=float)
+        eq.vacuum_toroidal_field.r0 = 1.0
+        eq.vacuum_toroidal_field.b0 = np.asarray([0.1, 0.2, 0.3], dtype=float)
         eq.time_slice.resize(3)
         for ts in eq.time_slice:
             ts.profiles_1d.psi = np.asarray([1.8, 1.9, 2.0, 2.1])
@@ -77,7 +82,7 @@ def interpolation_entry_path_directory(tmp_path_factory):
 
             ts.profiles_2d.resize(4)
             for p2d in ts.profiles_2d:
-                p2d.psi = np.asarray(np.random.rand(9, 3))
+                p2d.psi = np.asarray(rand_generator.random((9, 3)))
                 p2d.grid.dim1 = np.asarray([0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7], dtype=float)
                 p2d.grid.dim2 = np.asarray([1.0, 2.0, 3.0], dtype=float)
         entry.put(eq)
diff --git a/backend/tests/test_data_endpoints.py b/backend/tests/test_data_endpoints.py
index bad33101..ac5ba017 100644
--- a/backend/tests/test_data_endpoints.py
+++ b/backend/tests/test_data_endpoints.py
@@ -195,7 +195,7 @@ def test_plot_data_with_simple_operations_errors(entry_path):
 
 def test_plot_data_smoothing_with_wrong_target_node(entry_path):
     parameters = {
-        "uri": f"imas:hdf5?path={entry_path}#core_profiles/time",  # targetet quantity must be time-based
+        "uri": f"imas:hdf5?path={entry_path}#core_profiles/time",  # targeted quantity must be time-based
         "smoothing_method": "gaussian_filter",
         "gaussian_smoothing_sigma": 1,
     }
@@ -272,6 +272,57 @@ def test_plot_data_requires_gaussian_sigma(entry_path):
     assert "gaussian_smoothing_sigma is required" in response.text
 
 
+def test_combined_features(entry_path, interpolation_entry_path_directory):
+    """
+    Single test exercising all data manipulation features:
+    simple operations, smoothing, interpolation (exact_value), and signal operations.
+    """
+    # --- Part 1: simple ops + gaussian smoothing + signal ops (entry_path) ---
+    db = f"imas:hdf5?path={entry_path}"
+    parameters = {
+        "uri": f"{db}#core_profiles/global_quantities/ip",
+        "operations": ["add:2", "mul:3"],
+        "smoothing_method": "gaussian_filter",
+        "gaussian_smoothing_sigma": 1,
+        "signal_operations": [f"add:{db}#core_profiles/time"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+
+    response_body = response.json()
+    assert response_body["data"]["value"] == pytest.approx([11.28, 14.20, 18.0, 21.80, 24.72], 0.1)
+
+    # --- Part 2: different simple ops + savgol smoothing + exact_value interpolation + signal ops ---
+    db_names = [
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_1",
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_2",
+    ]
+
+    parameters = {
+        "uri": f"{db_names[0]}#equilibrium/vacuum_toroidal_field/b0",
+        "operations": ["mul:10", "pow:2"],
+        "smoothing_method": "savitzky-golay_filter",
+        "savgol_smoothing_window_length": 3,
+        "savgol_smoothing_polyorder": 1,
+        "signal_operations": [f"add:{db_names[1]}#equilibrium/vacuum_toroidal_field/b0"],
+        "interpolate_over": [f"{db_names[1]}#equilibrium/vacuum_toroidal_field/b0"],
+        "interpolation_method": "exact_value",
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+
+    response_body = response.json()
+    # db_1 b0: [0.1,0.2,0.3,0.4], mul:10+pow:2 -> [1,4,9,16]
+    # savgol wl=3 po=1 -> [0.6667,4.6667,9.6667,15.6667]
+    # exact_value interpolation on union [1,2,3,4] -> no change
+    # db_2 b0: [0.1,0.2,0.3]
+    # resampled to [1,2,3,4] with exact_value: [0.1,0.2,0.3,None]
+    # signal add propagates missing operand data: [0.7667,4.8667,9.9667,None]
+    values = response_body["data"]["value"]
+    assert values[:3] == pytest.approx([0.76, 4.86, 9.96], 0.01)
+    assert values[3] is None
+
+
 def test_plot_data_requires_savgol_window_length_and_polyorder(entry_path):
     base_parameters = {
         "uri": f"imas:hdf5?path={entry_path}#core_profiles/profiles_1d[:]/time",
@@ -320,3 +371,118 @@ def test_plot_data_coordinate_aliases(entry_path, expected_unit):
 
     assert dim2_coordinate["name"].lower() == "theta"
     assert dim2_coordinate["unit"].lower() == expected_unit[1]
+
+
+def test_plot_data_with_signal_operations_rejects_different_units(entry_path):
+    parameters = {
+        "uri": f"imas:hdf5?path={entry_path}#core_profiles/time",
+        "signal_operations": [f"add:imas:hdf5?path={entry_path}#core_profiles/global_quantities/ip"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 466
+    assert "Cannot add signals with different units" in response.json()["message"]
+
+
+@pytest.mark.parametrize(("operation", "expected_unit"), [("mul", "s*s"), ("div", "")])
+def test_plot_data_with_signal_operations_updates_unit(entry_path, operation, expected_unit):
+    uri = f"imas:hdf5?path={entry_path}#core_profiles/time"
+    response = pytest.test_client.get(
+        "/data/plot_data",
+        params={"uri": uri, "signal_operations": [f"{operation}:{uri}"]},
+    )
+
+    assert response.status_code == 200
+    assert response.json()["data"]["unit"] == expected_unit
+
+
+def test_plot_data_with_signal_operations_same_shape_different_uris(interpolation_entry_path_directory):
+    db_names = [
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_1",
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_2",
+    ]
+
+    parameters = {
+        "uri": f"{db_names[0]}#equilibrium/time_slice[0:2]/profiles_2d[0]/grid/dim2",
+        "signal_operations": [f"add:{db_names[1]}#equilibrium/time_slice[0:2]/profiles_2d[0]/grid/dim2"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+
+    response_body = response.json()
+    assert response_body["data"]["value"] == [[2.0, 4.0, 6.0], [2.0, 4.0, 6.0]]
+
+
+def test_plot_data_with_signal_operations_and_interpolation(interpolation_entry_path_directory):
+    db_names = [
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_1",
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_2",
+    ]
+
+    parameters = {
+        "uri": f"{db_names[0]}#equilibrium/time",
+        "signal_operations": [f"add:{db_names[1]}#equilibrium/time"],
+        "interpolate_over": [f"{db_names[1]}#equilibrium/time"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+
+    # time1: [1, 2, 3, 4]
+    # time2: [1, 2, 3]
+    response_body = response.json()
+    assert response_body["data"]["value"] == [2.0, 4.0, 6.0, None]
+    assert response_body["data"]["shape"] == [4]
+    assert response_body["data"]["downsampled_shape"] == [4]
+
+    # reversed order
+    parameters = {
+        "uri": f"{db_names[1]}#equilibrium/time",
+        "signal_operations": [f"add:{db_names[0]}#equilibrium/time"],
+        "interpolate_over": [f"{db_names[0]}#equilibrium/time"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+
+    # time1: [1, 2, 3]
+    # time2: [1, 2, 3, 4]
+    response_body = response.json()
+    assert response_body["data"]["value"] == [2.0, 4.0, 6.0, None]
+    # Interpolation expands the source from three to four samples.
+    assert response_body["data"]["shape"] == [4]
+    assert response_body["data"]["downsampled_shape"] == [4]
+
+
+def test_plot_data_with_signal_operations_and_interpolation_2d(interpolation_entry_path_directory):
+    db_names = [
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_1",
+        f"imas:hdf5?path={interpolation_entry_path_directory}/interpolation_db_2",
+    ]
+
+    parameters = {
+        "uri": f"{db_names[0]}#equilibrium/time_slice[:]/profiles_2d[:]/psi",
+        "signal_operations": [f"add:{db_names[1]}#equilibrium/time_slice[:]/profiles_2d[:]/psi"],
+        "interpolate_over": [f"{db_names[1]}#equilibrium/time_slice[:]/profiles_2d[:]/psi"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+    response_body = response.json()
+
+    data = np.array(response_body["data"]["value"], dtype=float)
+    # data shape reflects common coordinates (reversed): [time, profiles_2d, dim2, dim1]
+    assert data.shape == (4, 4, 3, 12)
+    # db_1 and db_2 have disjoint valid dim1 locations after interpolation, so operand NaNs propagate.
+    assert np.count_nonzero(~np.isnan(data)) == 0
+
+    # ---- reversed: db_2 primary, db_1 operand ----
+    parameters = {
+        "uri": f"{db_names[1]}#equilibrium/time_slice[:]/profiles_2d[:]/psi",
+        "signal_operations": [f"add:{db_names[0]}#equilibrium/time_slice[:]/profiles_2d[:]/psi"],
+        "interpolate_over": [f"{db_names[0]}#equilibrium/time_slice[:]/profiles_2d[:]/psi"],
+    }
+    response = pytest.test_client.get("/data/plot_data", params=parameters)
+    assert response.status_code == 200
+    response_body = response.json()
+
+    data = np.array(response_body["data"]["value"], dtype=float)
+    assert data.shape == (4, 4, 3, 12)
+    # db_1 and db_2 have disjoint valid dim1 locations after interpolation, so operand NaNs propagate.
+    assert np.count_nonzero(~np.isnan(data)) == 0
diff --git a/backend/tests/test_data_manipulation.py b/backend/tests/test_data_manipulation.py
index 5d637c2f..b4e47f37 100644
--- a/backend/tests/test_data_manipulation.py
+++ b/backend/tests/test_data_manipulation.py
@@ -1,12 +1,83 @@
 import numpy as np
 import pytest
+from ibex.data_source.exception import InvalidParametersException
 from ibex.data_source.imas_python_source_utils import (
     apply_gaussian_filter,
     apply_savgol_filter,
+    apply_signal_operations,
     apply_simple_operations,
+    combine_signal_units,
 )
 
 
+def test_apply_signal_operations_addition():
+    addend_uri = "imas:hdf5?path=/dummy/interpolation_db_1#equilibrium/time_slice[:]/profiles_2d[:]/psi"
+    data = np.array([1.0, 2.0, 3.0])
+    operations = [f"add:{addend_uri}"]
+    signal_data_by_uri = {addend_uri: np.array([10.0, 20.0, 30.0])}
+    result = apply_signal_operations(data, operations, signal_data_by_uri)
+    assert np.allclose(result, [11.0, 22.0, 33.0])
+
+
+def test_apply_signal_operations_subtraction():
+    subtrahend_uri = "imas:hdf5?path=/dummy/interpolation_db_2#equilibrium/time_slice[:]/profiles_2d[:]/psi"
+    data = np.array([10.0, 20.0, 30.0])
+    operations = [f"sub:{subtrahend_uri}"]
+    signal_data_by_uri = {subtrahend_uri: np.array([1.0, 2.0, 3.0])}
+    result = apply_signal_operations(data, operations, signal_data_by_uri)
+    assert np.allclose(result, [9.0, 18.0, 27.0])
+
+
+def test_apply_signal_operations_multiplication():
+    factor_uri = "imas:hdf5?path=/dummy/interpolation_db_1#equilibrium/time_slice[:]/profiles_2d[:]/psi"
+    data = np.array([1.0, 2.0, 3.0])
+    operations = [f"mul:{factor_uri}"]
+    signal_data_by_uri = {factor_uri: np.array([2.0, 3.0, 4.0])}
+    result = apply_signal_operations(data, operations, signal_data_by_uri)
+    assert np.allclose(result, [2.0, 6.0, 12.0])
+
+
+def test_apply_signal_operations_division():
+    divisor_uri = "imas:hdf5?path=/dummy/interpolation_db_2#equilibrium/time_slice[:]/profiles_2d[:]/psi"
+    data = np.array([10.0, 20.0, 30.0])
+    operations = [f"div:{divisor_uri}"]
+    signal_data_by_uri = {divisor_uri: np.array([2.0, 5.0, 6.0])}
+    result = apply_signal_operations(data, operations, signal_data_by_uri)
+    assert np.allclose(result, [5.0, 4.0, 5.0])
+
+
+def test_apply_signal_operations_preserves_operand_nans():
+    operand_uri = "imas:hdf5?path=/dummy/interpolation_db_1#equilibrium/time_slice[:]/profiles_2d[:]/psi"
+    data = np.array([10.0, 20.0, 30.0])
+    signal_data_by_uri = {operand_uri: np.array([1.0, np.nan, 3.0])}
+    result = apply_signal_operations(data, [f"add:{operand_uri}"], signal_data_by_uri)
+    assert np.allclose(result, [11.0, np.nan, 33.0], equal_nan=True)
+
+
+@pytest.mark.parametrize(
+    ("left_unit", "right_unit", "operation", "expected"),
+    [
+        ("kg", "kg", "add", "kg"),
+        ("kg", "kg", "sub", "kg"),
+        ("kg", "kg", "mul", "kg*kg"),
+        ("kg", "kg", "div", ""),
+        ("kg", "s", "div", "kg/s"),
+        ("", "kg", "mul", "kg"),
+    ],
+)
+def test_combine_signal_units(left_unit, right_unit, operation, expected):
+    assert combine_signal_units(left_unit, right_unit, operation) == expected
+
+
+@pytest.mark.parametrize("operation", ["add", "sub"])
+def test_combine_signal_units_rejects_different_units(operation):
+    with pytest.raises(
+        InvalidParametersException,
+        match=rf"Cannot {operation} signals with different units",
+    ):
+        combine_signal_units("kg", "s", operation)
+
+
 def test_apply_gaussian_smoothing():
     data = np.array([10.25, 12.8, 15.4, 18.15, 21.0, 24.35, 27.6, 30.2, 33.75, 36.1])
 
@@ -80,3 +151,31 @@ def test_apply_simple_operations_uses_order():
     # add then mul -> (5+1)*2 = 12
     result = apply_simple_operations(data, ["add:1", "mul:2"])
     assert result == pytest.approx([12.0])
+
+
+def test_apply_signal_operations_uses_order():
+    uri_a = "some/uri/a"
+    uri_b = "some/uri/b"
+    data = np.array([5.0])
+    signal_data_by_uri = {uri_a: np.array([2.0]), uri_b: np.array([1.0])}
+    # mul then add -> (5*2)+1 = 11
+    result = apply_signal_operations(data, [f"mul:{uri_a}", f"add:{uri_b}"], signal_data_by_uri)
+    assert result == pytest.approx([11.0])
+    # add then mul -> (5+1)*2 = 12
+    result = apply_signal_operations(data, [f"add:{uri_b}", f"mul:{uri_a}"], signal_data_by_uri)
+    assert result == pytest.approx([12.0])
+
+
+def test_apply_signal_operations_2D():
+    uri = "some/uri"
+    data = [np.array([1.0, 2.0]), np.array([3.0, 4.0])]
+    signal_data_by_uri = {uri: np.array([2.0, 3.0])}
+    result = apply_signal_operations(data, [f"add:{uri}", f"mul:{uri}"], signal_data_by_uri)
+    assert np.asarray(result[0]) == pytest.approx([6.0, 15.0])
+    assert np.asarray(result[1]) == pytest.approx([10.0, 21.0])
+
+
+def test_apply_signal_operations_rejects_division_by_zero():
+    uri = "some/uri"
+    with pytest.raises(InvalidParametersException, match="Division by zero is not allowed"):
+        apply_signal_operations(np.array([1.0]), [f"div:{uri}"], {uri: np.array([0.0])})
diff --git a/docs/source/developers_manual/backend_development/data_manipulation.rst b/docs/source/developers_manual/backend_development/data_manipulation.rst
index 39263c23..f8350d2b 100644
--- a/docs/source/developers_manual/backend_development/data_manipulation.rst
+++ b/docs/source/developers_manual/backend_development/data_manipulation.rst
@@ -14,10 +14,16 @@ The backend applies the manipulation stages in this order:
 1. simple data operations
 2. data smoothing
 3. data interpolation
-4. downsampling
+4. signal operations
+5. downsampling
 
 This means later stages operate on the output of earlier ones when the corresponding request parameters are enabled.
 
+.. note::
+
+   Signal operations are applied after data interpolation and before downsampling.
+   See the dedicated section below for details.
+
 Data smoothing
 ---------------
 
@@ -162,3 +168,69 @@ Order matters (multiply then add yields different result):
    curl -X 'GET' \
      '<IBEX_server_address>/data/plot_data?uri=<IMAS_URI>&operations=mul:3&operations=add:2' \
      -H 'accept: application/json'
+
+Signal operations
+------------------
+
+IBEX supports element-wise arithmetic between two signals (datasets) during a ``/data/plot_data`` request.
+Unlike simple scalar operations which use fixed numeric values, signal operations use **other IMAS signals** as operands.
+This enables combining data from different IMAS nodes arithmetically, for example subtracting a background signal from a measurement.
+
+Signal operations are applied **after data interpolation** and **before downsampling**.
+When a signal referenced in a signal operation was not already interpolated via the ``interpolate_over`` parameter, it is interpolated on the fly onto the common coordinate grid during this step.
+Operations are applied in the order they appear in the request.
+
+Configuration
+~~~~~~~~~~~~~~
+
+Signal operations are configured through the ``signal_operations`` parameter of the ``/data/plot_data/`` endpoint.
+It accepts a list of strings in the format ``type:uri``, where ``type`` is one of the supported operations and ``uri`` is the IMAS URI of the operand signal.
+
+The following operation types are supported:
+
+* ``add`` — addition
+* ``sub`` — subtraction
+* ``mul`` — multiplication
+* ``div`` — division
+* ``pow`` — exponentiation
+* ``root`` — Nth root
+
+Division by zero is rejected by the backend.
+
+The full list of available operations can be retrieved from the ``/info/data_manipulation_methods/`` endpoint.
+
+Implementation
+~~~~~~~~~~~~~~~
+
+Signal operations are applied in ``apply_signal_operations()`` in
+``backend/ibex/data_source/imas_python_source_utils.py``.
+
+The function iterates over the list of operation strings in order, looking up each
+operand signal's interpolated data from a dictionary populated during the interpolation step. For each operation the corresponding arithmetic function is applied element-wise between the current result and the operand signal data.
+
+The implementation handles both flat arrays and nested lists of arrays (higher-dimensional data) recursively.
+
+Example usage
+~~~~~~~~~~~~~~
+
+The following examples demonstrate how signal operations can be enabled for testing purposes.
+
+Addition of two signals:
+
+.. code-block:: bash
+
+   curl -X 'GET' \
+     '<IBEX_server_address>/data/plot_data?uri=<IMAS_URI>&signal_operations=add:<other_signal_URI>' \
+     -H 'accept: application/json'
+
+Multiple signal operations (subtraction then multiplication):
+
+.. code-block:: bash
+
+   curl -X 'GET' \
+     '<IBEX_server_address>/data/plot_data?uri=<IMAS_URI>&signal_operations=sub:<other_signal_URI>&signal_operations=mul:<yet_another_URI>' \
+     -H 'accept: application/json'
+
+
+
+