diff --git a/arithmetics-design/convention.md b/arithmetics-design/convention.md
index 64e3f16c..05c3f9c7 100644
--- a/arithmetics-design/convention.md
+++ b/arithmetics-design/convention.md
@@ -122,10 +122,10 @@ it does when no dimension matches. The same goes for a 4×4 array against
 `(a: 4, b: 4)`: sizes cannot tell `(a, b)` from `(b, a)`. To name the
 dimensions, wrap the array in a DataArray.
 
-> **TODO — not yet implemented ([#736]).** Today an unlabeled array pairs
-> with the *leading* dimensions positionally, which silently guesses in the
-> ambiguous cases above. The pairing rule builds on the `as_dataarray` /
-> coords-as-truth seam ([#732], merged — now `linopy.alignment`).
+A scalar broadcasts over every dimension and so needs no pairing. A 0-d
+array is treated as a scalar; a Python `list` is read as a numpy array
+(it carries values, not labels). Implemented in `linopy.alignment`
+([#736]).
 
 ### §8. Shared dimensions must match exactly
 
diff --git a/arithmetics-design/legacy-removal.md b/arithmetics-design/legacy-removal.md
index c19072be..5d563044 100644
--- a/arithmetics-design/legacy-removal.md
+++ b/arithmetics-design/legacy-removal.md
@@ -72,6 +72,11 @@ source tree.
   (`_project_onto_multiindex_levels`, `_LevelProjection`) stays:
   full-coverage full-level projections remain legal under v1 (they are
   the same coordinate spelled differently, §8).
+- `_dims_for_unlabeled_operand`: drop the legacy positional-pairing
+  fallback (the `warn_legacy(...)` branches plus the `return
+  list(candidates)`); the v1 size-pairing — the `is_v1()` block that
+  raises on ambiguity / no-match — becomes the whole function. The
+  `as_constant` / `_pair_axes_by_size` helpers stay (v1-clean).
 
 ### `linopy/piecewise.py` / `linopy/sos_reformulation.py`
 
diff --git a/linopy/alignment.py b/linopy/alignment.py
index 9649474c..919a5ae4 100644
--- a/linopy/alignment.py
+++ b/linopy/alignment.py
@@ -37,7 +37,31 @@
 from xarray.namedarray.utils import is_dict_like
 
 from linopy.constants import HELPER_DIMS
-from linopy.types import CoordsLike, DimsLike
+from linopy.types import UNLABELED_TYPES, CoordsLike, DimsLike
+
+
+def as_constant(other: Any) -> Any:
+    """
+    Normalize a degenerate operand for arithmetic on entry.
+
+    Two normalizations let the operators treat every numeric operand the
+    same way downstream:
+
+    - a Python ``list`` carries array data but no numeric operators
+      (``-[1, 2]`` is a ``TypeError``, ``[1, 2] * x`` repeats the list), so
+      it becomes a numpy array;
+    - a 0-d numpy array (``np.array(1)``) is unwrapped to a Python scalar so
+      it takes the scalar fast-path instead of size-pairing.
+
+    Everything else passes through unchanged — typed constants, DataArrays,
+    Variables, and Expressions all already behave.
+    """
+    if isinstance(other, list):
+        other = np.asarray(other)
+    if isinstance(other, np.ndarray) and other.ndim == 0:
+        return other.item()
+    return other
+
 
 if TYPE_CHECKING:
     from linopy.expressions import LinearExpression, QuadraticExpression
@@ -353,10 +377,13 @@ def as_dataarray(
         if isinstance(arr, np.number):
             arr = float(arr)
         if dims is None:
+            # A scalar broadcasts over the coords' dims, but never over a
+            # helper dim (e.g. ``_term``) — those are storage book-keeping,
+            # not user axes.
             if isinstance(coords, Coordinates):
-                dims = coords.dims
+                dims = [d for d in coords.dims if d not in HELPER_DIMS]
             elif is_dict_like(coords) and np.ndim(arr) == 0:
-                dims = list(coords.keys())
+                dims = [d for d in coords.keys() if d not in HELPER_DIMS]
         arr = DataArray(arr, coords=coords, dims=dims, **kwargs)
 
     elif not isinstance(arr, DataArray):
@@ -518,110 +545,223 @@ def _enforce_implicit_projections(projections: list[_LevelProjection]) -> None:
             )
 
 
-def _broadcast_to_coords(
-    arr: Any,
-    coords: CoordsLike | None = None,
-    dims: DimsLike | None = None,
-    **kwargs: Any,
-) -> tuple[DataArray, list[_LevelProjection]]:
+def _pair_axes_by_size(
+    shape: tuple[int, ...], sizes: dict[Hashable, int]
+) -> tuple[list[Hashable] | None, str | None]:
     """
-    Convert ``arr`` and broadcast it against ``coords`` (shared mechanics).
+    Pair each axis of an unlabeled array with the operand dim of matching size.
 
-    Returns the broadcast DataArray together with the MultiIndex-level
-    projections performed along the way, so the public entry points can
-    apply their own policy (warn or raise) to partial projections and
-    coverage gaps.
+    The pairing must be determined by the sizes alone (v1 convention,
+    coordinate-alignment intro): every axis size must match exactly one
+    operand dim, and no two axes may share a size. Returns
+    ``(dims, None)`` on success or ``(None, problem)`` where ``problem``
+    describes why the pairing is impossible or ambiguous.
     """
-    if coords is None:
-        return as_dataarray(arr, coords, dims, **kwargs), []
+    by_size: dict[int, list[Hashable]] = {}
+    for d, n in sizes.items():
+        by_size.setdefault(n, []).append(d)
+
+    axes_per_size: dict[int, int] = {}
+    for s in shape:
+        axes_per_size[s] = axes_per_size.get(s, 0) + 1
+
+    for s, n_axes in axes_per_size.items():
+        candidates = by_size.get(s, [])
+        if len(candidates) < n_axes:
+            return None, (
+                f"no unambiguous dimension match for an axis of length {s}: "
+                f"the operand has dimensions {dict(sizes)}."
+            )
+        if len(candidates) > 1 or n_axes > 1:
+            return None, (
+                f"axis of length {s} could pair with any of "
+                f"{sorted(candidates, key=str)} — sizes alone cannot decide."
+            )
 
-    if isinstance(coords, list | tuple) and any(isinstance(c, tuple) for c in coords):
-        # xarray reads bare `(a, b)` as `(dim_name, values)`; normalize so a
-        # coords entry passed as a tuple behaves identically to a list.
-        coords = [list(c) if isinstance(c, tuple) else c for c in coords]
+    return [by_size[s][0] for s in shape], None
 
-    expected = _coords_to_dict(coords, dims=dims)
-    if not expected:
-        return as_dataarray(arr, coords, dims, **kwargs), []
 
-    if isinstance(arr, pd.Series | pd.DataFrame):
-        converted = _named_pandas_to_dataarray(arr)
-        if converted is not None:
-            arr = converted
+def _dims_for_unlabeled_operand(
+    shape: tuple[int, ...], expected: dict[Hashable, Any]
+) -> list[Hashable]:
+    """
+    Choose dim names for an unlabeled (numpy / list / polars) input.
 
-    if not isinstance(arr, DataArray):
-        # numpy/polars/unnamed-pandas inputs are positional — their only
-        # meaningful information is the values; any axis labels are
-        # auto-generated. Default dims to coords' keys so the conversion
-        # labels axes correctly (instead of dim_0/dim_1), then re-assign
-        # coords from expected so positional inputs align to coords by
-        # position. A shape mismatch surfaces here as a clear xarray
-        # "conflicting sizes" error rather than a confusing
-        # "coordinates do not match" further down.
-        if dims is None:
-            dims = list(expected)
-        arr = as_dataarray(arr, coords, dims=dims, **kwargs)
-        # Skip MultiIndex dims — re-assigning a PandasMultiIndex coord emits
-        # a FutureWarning and isn't needed (the conversion already used it).
-        arr = arr.assign_coords(
-            {
-                d: expected[d]
-                for d in arr.dims
-                if d in expected and not isinstance(arr.indexes.get(d), pd.MultiIndex)
-            }
+    Used everywhere an unlabeled array meets a known set of dims — bounds
+    and masks in ``add_variables`` / ``add_constraints``, and arithmetic
+    operands (#736).
+
+    v1 (convention, coordinate-alignment intro): axes pair with the dims by
+    size; ambiguity or a missing match raises, with wrap-in-a-DataArray as
+    the documented resolution. Legacy: axes pair with the leading dims
+    positionally; a deprecation warning fires whenever the v1 pairing would
+    differ from or reject the positional one.
+    """
+    from linopy.semantics import is_v1, warn_legacy
+
+    # A 0-d operand has no axes to pair — it broadcasts over every dim, so it
+    # carries no dim names (matching a bare scalar).
+    if len(shape) == 0:
+        return []
+
+    # Helper dims (e.g. ``_term``) are storage book-keeping, never user axes,
+    # so they are not pairing candidates.
+    candidates = {d: v for d, v in expected.items() if d not in HELPER_DIMS}
+    sizes = {d: len(_as_index(v)) for d, v in candidates.items()}
+    paired, problem = _pair_axes_by_size(shape, sizes)
+    positional = list(candidates)[: len(shape)]
+
+    if is_v1():
+        if problem is not None:
+            raise ValueError(
+                f"Cannot pair an unlabeled array of shape {tuple(shape)} with "
+                f"the operand's dimensions: {problem} Wrap the array in an "
+                f"xarray.DataArray with explicit dims to name its axes."
+            )
+        assert paired is not None
+        return paired
+
+    # LEGACY: remove at 1.0 — positional pairing plus the transition warning.
+    if problem is not None:
+        warn_legacy(
+            f"An unlabeled array of shape {tuple(shape)} was paired with the "
+            f"operand's leading dimension(s) {positional} by position. Under "
+            f"the v1 convention this raises: {problem} Wrap the array in an "
+            f"xarray.DataArray with explicit dims to keep it working."
+        )
+    elif paired != positional:
+        warn_legacy(
+            f"An unlabeled array of shape {tuple(shape)} was paired with the "
+            f"operand's leading dimension(s) {positional} by position. Under "
+            f"the v1 convention it pairs by size instead — with {paired} — "
+            f"which gives a different result. Wrap the array in an "
+            f"xarray.DataArray with explicit dims to make the pairing explicit."
         )
+    return positional
 
-    arr, projections = _project_onto_multiindex_levels(arr, expected)
 
+def _matmul_operand_to_dataarray(
+    other: Any, coords: Coordinates, coord_dims: tuple[Hashable, ...]
+) -> DataArray:
+    """
+    Convert a non-expression ``@`` operand, pairing unlabeled axes by size.
+
+    Shared by ``LinearExpression.__matmul__`` and
+    ``QuadraticExpression.__matmul__``: :func:`_dims_for_positional_input`
+    decides which dims the contraction collapses (#736) — by size under v1,
+    positionally with a warning under legacy. Unlike the broadcast pipeline
+    this only converts (no reindex / expand / transpose).
+    """
+    expected = {d: coords[d] for d in coord_dims}
+    dims = _dims_for_positional_input(other, expected, None)
+    return as_dataarray(other, coords=coords, dims=dims)
+
+
+def _dims_for_positional_input(
+    arr: Any, expected: dict[Hashable, Any], dims: DimsLike | None
+) -> DimsLike | None:
+    """
+    Resolve the dim names a non-DataArray input's axes adopt.
+
+    An explicit ``dims`` is honored as given. Otherwise an unlabeled
+    array (numpy / list / polars) pairs its axes with ``expected`` by
+    size (#736); any other input falls back to the coords dims, minus
+    helper dims like ``_term`` which are never user axes.
+    """
+    if dims is not None:
+        return dims
+    if isinstance(arr, UNLABELED_TYPES) and np.ndim(arr) >= 1:
+        return _dims_for_unlabeled_operand(np.shape(arr), expected)
+    return [d for d in expected if d not in HELPER_DIMS]
+
+
+def _label_input(
+    arr: Any, expected: dict[Hashable, Any], dims: DimsLike | None, **kwargs: Any
+) -> DataArray:
+    """
+    Convert a non-DataArray input to a DataArray labelled against ``expected``.
+
+    The converter is handed ``expected`` (the normalized name→values dict),
+    not the raw sequence-form coords, so it selects coords by name — a
+    sequence would zip dims to coords by position, which is wrong once
+    size-pairing has chosen a non-leading dim.
+    """
+    dims = _dims_for_positional_input(arr, expected, dims)
+    arr = as_dataarray(arr, expected, dims=dims, **kwargs)
+    # Re-assign non-MultiIndex coords from ``expected`` (a MultiIndex coord
+    # re-assignment emits a FutureWarning and the conversion already used it).
+    return arr.assign_coords(
+        {
+            d: expected[d]
+            for d in arr.dims
+            if d in expected and not isinstance(arr.indexes.get(d), pd.MultiIndex)
+        }
+    )
+
+
+def _reindex_reordered_dims(arr: DataArray, expected: dict[Hashable, Any]) -> DataArray:
+    """
+    Reindex shared dims whose values match ``expected`` in a different order.
+
+    Disagreeing value *sets* are left for downstream xarray alignment;
+    only a pure reordering of the same values is conformed here.
+    """
     for dim, coord_values in expected.items():
-        if dim not in arr.dims:
-            continue
-        if isinstance(arr.indexes.get(dim), pd.MultiIndex):
+        if dim not in arr.dims or isinstance(arr.indexes.get(dim), pd.MultiIndex):
             continue
         expected_idx = _as_index(coord_values)
         actual_idx = arr.coords[dim].to_index()
         if actual_idx.equals(expected_idx):
             continue
-        # Same values, different order → reindex to match expected order.
-        # Different value sets are left alone for downstream xarray alignment.
         if len(actual_idx) == len(expected_idx) and set(actual_idx) == set(
             expected_idx
         ):
             arr = arr.reindex({dim: expected_idx})
+    return arr
+
 
-    # expand_dims prepends new dimensions and their coordinate variables;
-    # the subsequent transpose restores coords order. Both are no-ops when
-    # the array already matches. Reconstruct so the DataArray's coords
-    # iteration order also follows coords (a Dataset built from this picks
-    # up its dim order from coord insertion).
+def _expand_missing_dims(arr: DataArray, expected: dict[Hashable, Any]) -> DataArray:
+    """
+    Broadcast ``arr`` over ``expected`` dims it does not yet carry.
+
+    A MultiIndex-backed dim is broadcast against a proper ``Coordinates``
+    template: plain ``expand_dims`` would drop its level coords and leave a
+    degenerate flat index that fails to align downstream. The exception is
+    when ``arr`` already carries one of the MultiIndex's level names —
+    broadcasting would then raise on the conflicting index, so fall back to
+    ``expand_dims``.
+    """
     expand = {k: v for k, v in expected.items() if k not in arr.dims}
-    if expand:
-        # expand_dims drops the level coords of a MultiIndex-backed dim,
-        # leaving a degenerate flat index that fails to align downstream.
-        # Broadcast against a proper Coordinates template instead.
-        plain = {}
-        for dim, coord_values in expand.items():
-            mi = _as_multiindex(coord_values)
-            # Fall back to expand_dims when arr already carries one of the
-            # MultiIndex's level names as its own coord: broadcasting against
-            # the level coords would raise on the conflicting index.
-            if mi is None or set(mi.names) & (set(arr.coords) | set(arr.dims)):
-                plain[dim] = coord_values
-                continue
-            template = DataArray(
-                np.zeros(len(mi)),
-                coords=Coordinates.from_pandas_multiindex(mi, dim),
-                dims=[dim],
-            )
-            arr, _ = broadcast(arr, template)
-        if plain:
-            arr = arr.expand_dims(plain)
+    if not expand:
+        return arr
+    plain = {}
+    for dim, coord_values in expand.items():
+        mi = _as_multiindex(coord_values)
+        if mi is None or set(mi.names) & (set(arr.coords) | set(arr.dims)):
+            plain[dim] = coord_values
+            continue
+        template = DataArray(
+            np.zeros(len(mi)),
+            coords=Coordinates.from_pandas_multiindex(mi, dim),
+            dims=[dim],
+        )
+        arr, _ = broadcast(arr, template)
+    if plain:
+        arr = arr.expand_dims(plain)
+    return arr
+
+
+def _order_like_coords(arr: DataArray, expected: dict[Hashable, Any]) -> DataArray:
+    """
+    Transpose ``arr`` to ``coords`` dim order, then match coord iteration order.
 
+    The reconstruction makes a Dataset built from ``arr`` pick up its dim
+    order from coord insertion, not just the transpose.
+    """
     target_dims = tuple(d for d in expected if d in arr.dims) + tuple(
         d for d in arr.dims if d not in expected
     )
     arr = arr.transpose(*target_dims)
-
     coord_order = [c for c in target_dims if c in arr.coords] + [
         c for c in arr.coords if c not in target_dims
     ]
@@ -631,7 +771,48 @@ def _broadcast_to_coords(
             coords={c: arr.coords[c] for c in coord_order},
             name=arr.name,
         )
+    return arr
+
 
+def _broadcast_to_coords(
+    arr: Any,
+    coords: CoordsLike | None = None,
+    dims: DimsLike | None = None,
+    **kwargs: Any,
+) -> tuple[DataArray, list[_LevelProjection]]:
+    """
+    Convert ``arr`` and broadcast it against ``coords`` (shared mechanics).
+
+    Returns the broadcast DataArray together with the MultiIndex-level
+    projections performed along the way, so the public entry points can
+    apply their own policy (warn or raise) to partial projections and
+    coverage gaps. Unlabeled inputs pair their axes with the coords dims by
+    size (#736); see :func:`_label_input`.
+    """
+    if coords is None:
+        return as_dataarray(arr, coords, dims, **kwargs), []
+
+    if isinstance(coords, list | tuple) and any(isinstance(c, tuple) for c in coords):
+        # xarray reads bare `(a, b)` as `(dim_name, values)`; normalize so a
+        # tuple coords entry behaves identically to a list.
+        coords = [list(c) if isinstance(c, tuple) else c for c in coords]
+
+    expected = _coords_to_dict(coords, dims=dims)
+    if not expected:
+        return as_dataarray(arr, coords, dims, **kwargs), []
+
+    if isinstance(arr, pd.Series | pd.DataFrame):
+        converted = _named_pandas_to_dataarray(arr)
+        if converted is not None:
+            arr = converted
+
+    if not isinstance(arr, DataArray):
+        arr = _label_input(arr, expected, dims, **kwargs)
+
+    arr, projections = _project_onto_multiindex_levels(arr, expected)
+    arr = _reindex_reordered_dims(arr, expected)
+    arr = _expand_missing_dims(arr, expected)
+    arr = _order_like_coords(arr, expected)
     return arr, projections
 
 
diff --git a/linopy/examples.py b/linopy/examples.py
index 6e1cfb15..549574bf 100644
--- a/linopy/examples.py
+++ b/linopy/examples.py
@@ -2,6 +2,7 @@
 This module contains examples of linear programming models using the linopy library.
 """
 
+import xarray as xr
 from numpy import arange
 
 from linopy import Model
@@ -73,7 +74,7 @@ def benchmark_model(n: int = 10, integerlabels: bool = False) -> Model:
         naxis, maxis = [arange(n, dtype=float), arange(n).astype(str)]
     x = m.add_variables(coords=[naxis, maxis])
     y = m.add_variables(coords=[naxis, maxis])
-    m.add_constraints(x - y >= naxis)
+    m.add_constraints(x - y >= xr.DataArray(naxis, dims=["dim_0"]))
     m.add_constraints(x + y >= 0)
     m.add_objective((2 * x).sum() + y.sum())
     return m
diff --git a/linopy/expressions.py b/linopy/expressions.py
index b7a166d1..ed526a1c 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -44,7 +44,13 @@
 from types import EllipsisType, NotImplementedType
 
 from linopy import constraints, variables
-from linopy.alignment import as_dataarray, broadcast_to_coords, fill_missing_coords
+from linopy.alignment import (
+    _matmul_operand_to_dataarray,
+    as_constant,
+    as_dataarray,
+    broadcast_to_coords,
+    fill_missing_coords,
+)
 from linopy.common import (
     EmptyDeprecationWrapper,
     LocIndexer,
@@ -680,9 +686,7 @@ def _add_constant_v1(
             if isinstance(other, float) and np.isnan(other):
                 check_user_nan()
             return self.assign(const=self.const + other)
-        da = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        da = broadcast_to_coords(other, coords=self.coords, strict=False)
         if da.isnull().any():
             check_user_nan()
         self_const, da, needs_data_reindex = self._align_constant(
@@ -708,9 +712,7 @@ def _add_constant_legacy(
             if isinstance(other, float) and np.isnan(other):
                 check_user_nan()
             return self.assign(const=self.const.fillna(0) + other)
-        da = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        da = broadcast_to_coords(other, coords=self.coords, strict=False)
         if da.isnull().any():
             check_user_nan()
         self_const, da, needs_data_reindex = self._align_constant(
@@ -752,9 +754,7 @@ def _apply_constant_op_v1(
         # §5: user NaN raised before we get here.
         if isinstance(other, float) and np.isnan(other):
             check_user_nan(op_kind=op_kind)
-        factor = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        factor = broadcast_to_coords(other, coords=self.coords, strict=False)
         if factor.isnull().any():
             check_user_nan(op_kind=op_kind)
         self_const, factor, needs_data_reindex = self._align_constant(
@@ -785,9 +785,7 @@ def _apply_constant_op_legacy(
         # factor → fill_value (0 for mul, 1 for div), coeffs/const → 0.
         if isinstance(other, float) and np.isnan(other):
             check_user_nan(op_kind=op_kind)
-        factor = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        factor = broadcast_to_coords(other, coords=self.coords, strict=False)
         if factor.isnull().any():
             check_user_nan(op_kind=op_kind)
         self_const, factor, needs_data_reindex = self._align_constant(
@@ -822,6 +820,7 @@ def _divide_by_constant(
         )
 
     def __div__(self: GenericExpression, other: SideLike) -> GenericExpression:
+        other = as_constant(other)
         try:
             if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
                 raise TypeError(
@@ -1280,6 +1279,7 @@ def to_constraint(
         which are moved to the left-hand-side and constant values which are moved
         to the right-hand side.
         """
+        rhs = as_constant(rhs)
         if self.is_constant and is_constant(rhs):
             raise ValueError(
                 f"Both sides of the constraint are constant. At least one side must contain variables. {self} {rhs}"
@@ -1294,9 +1294,7 @@ def to_constraint(
             # through ``sub`` into the RHS and reaches downstream
             # auto-mask handling as "no constraint at this row" (§12).
             if isinstance(rhs, CONSTANT_TYPES):
-                rhs = broadcast_to_coords(
-                    rhs, coords=self.coords, dims=self.coord_dims, strict=False
-                )
+                rhs = broadcast_to_coords(rhs, coords=self.coords, strict=False)
                 extra_dims = set(rhs.dims) - set(self.coord_dims)
                 if extra_dims:
                     logger.warning(
@@ -1320,9 +1318,7 @@ def to_constraint(
         # part of normal arithmetic, so we restore the original NaN mask
         # afterward).
         if isinstance(rhs, CONSTANT_TYPES):
-            rhs = broadcast_to_coords(
-                rhs, coords=self.coords, dims=self.coord_dims, strict=False
-            )
+            rhs = broadcast_to_coords(rhs, coords=self.coords, strict=False)
 
             extra_dims = set(rhs.dims) - set(self.coord_dims)
             if extra_dims:
@@ -1844,6 +1840,7 @@ def __add__(
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         if isinstance(other, QuadraticExpression):
             return other.__add__(self)
 
@@ -1879,6 +1876,7 @@ def __sub__(
         | LinearExpression
         | QuadraticExpression,
     ) -> LinearExpression | QuadraticExpression:
+        other = as_constant(other)
         try:
             return self.__add__(-other)
         except TypeError:
@@ -1903,6 +1901,7 @@ def __mul__(
         """
         Multiply the expr by a factor.
         """
+        other = as_constant(other)
         if isinstance(other, QuadraticExpression):
             return other.__rmul__(self)
 
@@ -1948,8 +1947,9 @@ def __matmul__(
         """
         Matrix multiplication with other, similar to xarray dot.
         """
+        other = as_constant(other)
         if not isinstance(other, LinearExpression | variables.Variable):
-            other = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+            other = _matmul_operand_to_dataarray(other, self.coords, self.coord_dims)
 
         common_dims = list(set(self.coord_dims).intersection(other.dims))
         return (self * other).sum(dim=common_dims)
@@ -2360,6 +2360,7 @@ def __mul__(self, other: SideLike) -> QuadraticExpression:
         """
         Multiply the expr by a factor.
         """
+        other = as_constant(other)
         if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
             raise TypeError(
                 "unsupported operand type(s) for *: "
@@ -2381,6 +2382,7 @@ def __add__(self, other: SideLike) -> QuadraticExpression:
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         try:
             if isinstance(other, CONSTANT_TYPES):
                 return self._add_constant(other)
@@ -2407,6 +2409,7 @@ def __sub__(self, other: SideLike) -> QuadraticExpression:
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         try:
             return self.__add__(-other)
         except TypeError:
@@ -2430,12 +2433,13 @@ def __matmul__(
         """
         Matrix multiplication with other, similar to xarray dot.
         """
+        other = as_constant(other)
         if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
             raise TypeError(
                 "Higher order non-linear expressions are not yet supported."
             )
 
-        other = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+        other = _matmul_operand_to_dataarray(other, self.coords, self.coord_dims)
         common_dims = list(set(self.coord_dims).intersection(other.dims))
         return (self * other).sum(dim=common_dims)
 
diff --git a/linopy/types.py b/linopy/types.py
index 6b4cf712..21a8dddd 100644
--- a/linopy/types.py
+++ b/linopy/types.py
@@ -39,6 +39,10 @@
     | pl.Series
 )
 CONSTANT_TYPES: tuple[type, ...] = get_args(ConstantLike)
+
+UnlabeledLike: TypeAlias = numpy.ndarray | list | pl.Series
+UNLABELED_TYPES = get_args(UnlabeledLike)
+
 SignLike: TypeAlias = str | numpy.ndarray | DataArray | Series | DataFrame
 MaskLike: TypeAlias = numpy.ndarray | DataArray | Series | DataFrame
 PathLike: TypeAlias = str | Path
diff --git a/linopy/variables.py b/linopy/variables.py
index bc4b1e8b..8c83a305 100644
--- a/linopy/variables.py
+++ b/linopy/variables.py
@@ -333,9 +333,7 @@ def to_linexpr(
         linopy.LinearExpression
             Linear expression with the variables and coefficients.
         """
-        coefficient = broadcast_to_coords(
-            coefficient, coords=self.coords, dims=self.dims, strict=False
-        )
+        coefficient = broadcast_to_coords(coefficient, coords=self.coords, strict=False)
         # §5: user-supplied NaN in the coefficient must raise (v1) / warn
         # (legacy) — it's the multiplicative analogue of ``x + nan_data``
         # and otherwise enters the expression silently. The default
diff --git a/test/test_alignment.py b/test/test_alignment.py
index d571892c..118b0d20 100644
--- a/test/test_alignment.py
+++ b/test/test_alignment.py
@@ -600,6 +600,25 @@ def test_extra_coord_entries_broadcast_in(self) -> None:
         assert list(da.coords["dim_0"].values) == ["a", "b"]
         assert list(da.coords["dim_2"].values) == ["A", "B"]
 
+    @pytest.mark.v1
+    def test_explicit_dims_bypass_size_pairing(self) -> None:
+        """
+        An explicit ``dims`` is honored positionally, like xarray — even when
+        size-pairing would otherwise be ambiguous (both coords dims size 4).
+        Pins the "infer order only when the user didn't name it" rule.
+        """
+        coords = {"a": [0, 1, 2, 3], "b": [4, 5, 6, 7]}  # both size 4
+
+        # No dims → size-pairing can't decide → raises.
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            broadcast_to_coords(np.arange(4), coords=coords, strict=False)
+
+        # Explicit dims=['a'] → the axis is labeled 'a' positionally, no
+        # pairing, no raise (matches xarray's positional dims assignment).
+        da = broadcast_to_coords(np.arange(4), coords=coords, dims=["a"], strict=False)
+        assert set(da.dims) == {"a", "b"}
+        assert (da.sel(b=4).values == np.arange(4)).all()
+
 
 # ---------------------------------------------------------------------------
 # Implicit MultiIndex-level projection — the legacy/v1 fork point
diff --git a/test/test_constraints.py b/test/test_constraints.py
index 27152f98..2de429d1 100644
--- a/test/test_constraints.py
+++ b/test/test_constraints.py
@@ -144,7 +144,6 @@ def test_constraint_assignment_with_reindex() -> None:
 @pytest.mark.parametrize(
     "rhs_factory",
     [
-        pytest.param(lambda m, v: v, id="numpy"),
         pytest.param(lambda m, v: xr.DataArray(v, dims=["dim_0"]), id="dataarray"),
         pytest.param(lambda m, v: pd.Series(v, index=v), id="series"),
         pytest.param(
@@ -168,10 +167,37 @@ def test_constraint_rhs_lower_dim(rhs_factory: Any) -> None:
     assert c.shape == (10, 10)
 
 
+@pytest.mark.legacy
+def test_constraint_rhs_unlabeled_lower_dim_legacy() -> None:
+    # An unlabeled array rhs pairs positionally with the leading dim under
+    # legacy; both dims are size 10 so the pairing is a size-coincidence.
+    m = Model()
+    naxis = np.arange(10, dtype=float)
+    maxis = np.arange(10).astype(str)
+    x = m.add_variables(coords=[naxis, maxis])
+    y = m.add_variables(coords=[naxis, maxis])
+
+    c = m.add_constraints(x - y >= naxis)
+    assert c.shape == (10, 10)
+
+
+@pytest.mark.v1
+def test_constraint_rhs_unlabeled_lower_dim_ambiguous_raises_v1() -> None:
+    # v1: both dims are size 10, so an unlabeled length-10 rhs cannot be
+    # paired by size — it raises (wrap in a DataArray to disambiguate).
+    m = Model()
+    naxis = np.arange(10, dtype=float)
+    maxis = np.arange(10).astype(str)
+    x = m.add_variables(coords=[naxis, maxis])
+    y = m.add_variables(coords=[naxis, maxis])
+
+    with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+        m.add_constraints(x - y >= naxis)
+
+
 @pytest.mark.parametrize(
     "rhs_factory",
     [
-        pytest.param(lambda m: np.ones((5, 3)), id="numpy"),
         pytest.param(lambda m: pd.DataFrame(np.ones((5, 3))), id="dataframe"),
     ],
 )
@@ -186,6 +212,29 @@ def test_constraint_rhs_higher_dim_constant_warns(
     assert "dimensions" in caplog.text
 
 
+@pytest.mark.legacy
+def test_constraint_rhs_unlabeled_higher_dim_warns_legacy(caplog: Any) -> None:
+    # Legacy: an unlabeled (5, 3) rhs pairs axis 0 with dim_0 positionally and
+    # broadcasts the extra axis, warning about the extra dimension.
+    m = Model()
+    x = m.add_variables(coords=[range(5)], name="x")
+
+    with caplog.at_level("WARNING", logger="linopy.expressions"):
+        m.add_constraints(x >= np.ones((5, 3)))
+    assert "dimensions" in caplog.text
+
+
+@pytest.mark.v1
+def test_constraint_rhs_unlabeled_higher_dim_raises_v1() -> None:
+    # v1: the (5, 3) rhs has an axis of length 3 matching no dim of x — it
+    # cannot be paired by size and raises.
+    m = Model()
+    x = m.add_variables(coords=[range(5)], name="x")
+
+    with pytest.raises(ValueError, match=r"no unambiguous dimension match"):
+        m.add_constraints(x >= np.ones((5, 3)))
+
+
 def test_constraint_rhs_higher_dim_dataarray_reindexes() -> None:
     """DataArray RHS with extra dims reindexes to expression coords (no raise)."""
     m = Model()
diff --git a/test/test_legacy_violations.py b/test/test_legacy_violations.py
index 1ccb1714..5482770a 100644
--- a/test/test_legacy_violations.py
+++ b/test/test_legacy_violations.py
@@ -34,12 +34,17 @@ class corresponds to a section of ``arithmetics-design/convention.md``
     §11 Non-dim coord conflict raises (v1)             → #295
     §11 Non-conflicting aux coords propagate through arithmetic
 
+Slice H — unlabeled-operand pairing (coordinate-alignment intro):
+    Unlabeled operands (numpy / list / polars) pair with the linopy
+    operand's dims by size; ambiguity or no-match raises (v1)   → #736
+
 Slice H — object scope (convention preamble):
     Non-linopy operands behave exactly like constant-only expressions
 """
 
 from __future__ import annotations
 
+import contextlib
 import operator
 import warnings
 from collections.abc import Generator
@@ -47,6 +52,7 @@ class corresponds to a section of ``arithmetics-design/convention.md``
 
 import numpy as np
 import pandas as pd
+import polars as pl
 import pytest
 import xarray as xr
 
@@ -167,6 +173,226 @@ def test_mul_broadcast_introduces_new_dim(self, x: Variable) -> None:
         assert set(result.coeffs.dims) == {"time", "scenario", "_term"}
 
 
+# =====================================================================
+# Coordinate-alignment intro — unlabeled operands pair by size (#736)
+# =====================================================================
+
+
+# The three unlabeled types must behave identically; each constructor takes a
+# value sequence so the same cases parametrize over numpy / list / polars.
+UNLABELED_1D = [
+    pytest.param(lambda v: np.asarray(v, dtype=float), id="numpy"),
+    pytest.param(lambda v: [float(x) for x in v], id="list"),
+    pytest.param(lambda v: pl.Series([float(x) for x in v]), id="polars"),
+]
+
+
+class TestUnlabeledPairing:
+    """
+    Unlabeled operands (numpy arrays, lists, polars Series) carry no labels,
+    so they pair with the linopy operand's dims *by size*. Ambiguity (two
+    dims share the size, or the array is square) or no size match raises
+    under v1; legacy pairs with the leading dims positionally and warns when
+    the v1 pairing would differ or reject.
+    """
+
+    @pytest.fixture
+    def xy(self) -> Variable:
+        # dims of distinct sizes so a 1-d operand pairs unambiguously
+        m = Model()
+        return m.add_variables(
+            coords=[pd.RangeIndex(3, name="a"), pd.RangeIndex(4, name="b")], name="xy"
+        )
+
+    @pytest.fixture
+    def square(self) -> Variable:
+        # both dims size 4 → a 1-d length-4 operand is ambiguous
+        m = Model()
+        return m.add_variables(
+            coords=[pd.RangeIndex(4, name="p"), pd.RangeIndex(4, name="q")], name="sq"
+        )
+
+    @pytest.fixture
+    def wide(self) -> Variable:
+        # four dims of distinct sizes — a lower-rank operand pairs a subset
+        m = Model()
+        return m.add_variables(
+            coords=[
+                pd.RangeIndex(3, name="a"),
+                pd.RangeIndex(4, name="b"),
+                pd.RangeIndex(5, name="c"),
+                pd.RangeIndex(6, name="d"),
+            ],
+            name="wide",
+        )
+
+    # -- 1-d operands -----------------------------------------------------
+
+    @pytest.mark.v1
+    @pytest.mark.parametrize("make", UNLABELED_1D)
+    def test_v1_pairs_by_size(self, xy: Variable, make: Any) -> None:
+        # length-4 array pairs with dim "b" (size 4), not the leading "a" (3)
+        result = (1 * xy) + make(range(4))
+        assert set(result.const.dims) == {"a", "b"}
+        assert result.const.sizes == {"a": 3, "b": 4}
+
+    @pytest.mark.v1
+    @pytest.mark.parametrize("make", UNLABELED_1D)
+    def test_v1_ambiguous_square_raises(self, square: Variable, make: Any) -> None:
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            (1 * square) + make(range(4))
+
+    @pytest.mark.v1
+    @pytest.mark.parametrize("make", UNLABELED_1D)
+    def test_v1_no_size_match_raises(self, xy: Variable, make: Any) -> None:
+        with pytest.raises(ValueError, match=r"no unambiguous dimension match"):
+            (1 * xy) + make(range(7))
+
+    @pytest.mark.v1
+    def test_v1_dataarray_wrapping_resolves_ambiguity(self, square: Variable) -> None:
+        # the documented escape hatch: name the axis with a DataArray
+        result = (1 * square) + xr.DataArray(np.arange(4.0), dims=["p"])
+        assert set(result.const.dims) == {"p", "q"}
+
+    # -- multi-dim operands (numpy only — list / polars are 1-d) ----------
+
+    @pytest.mark.v1
+    def test_v1_size_order_independent(self, xy: Variable) -> None:
+        # a 2-d (4, 3) operand pairs (b, a) by size regardless of axis order.
+        result = (1 * xy) + np.ones((4, 3))
+        assert result.const.sizes == {"a": 3, "b": 4}
+
+    @pytest.mark.v1
+    def test_v1_multidim_square_ambiguous_raises(self, square: Variable) -> None:
+        # a 2-d (4, 4) operand against dims (p: 4, q: 4) — sizes cannot tell
+        # (p, q) from (q, p). Exercises the multi-axis-same-size branch the 1-d
+        # cases never reach.
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            (1 * square) + np.ones((4, 4))
+
+    @pytest.mark.v1
+    def test_v1_multidim_no_size_match_raises(self, xy: Variable) -> None:
+        # a 2-d (5, 3) operand against (a: 3, b: 4) — the length-5 axis matches
+        # no dim.
+        with pytest.raises(ValueError, match=r"no unambiguous dimension match"):
+            (1 * xy) + np.ones((5, 3))
+
+    @pytest.mark.v1
+    def test_v1_lower_rank_operand_pairs_subset_and_broadcasts(
+        self, wide: Variable
+    ) -> None:
+        # a 2-d (4, 5) operand against four dims pairs its axes with (b, c) by
+        # size and broadcasts over the unpaired (a, d).
+        result = (1 * wide) + np.ones((4, 5))
+        assert result.const.sizes == {"a": 3, "b": 4, "c": 5, "d": 6}
+
+    @pytest.mark.v1
+    def test_v1_ambiguous_axis_within_higher_rank_raises(self) -> None:
+        # a 2-d (4, 5) operand where the length-4 axis matches two dims of the
+        # operand (a: 4, b: 4) — ambiguous even though the length-5 axis is
+        # unique.
+        m = Model()
+        y = m.add_variables(
+            coords=[
+                pd.RangeIndex(4, name="a"),
+                pd.RangeIndex(4, name="b"),
+                pd.RangeIndex(5, name="c"),
+            ],
+            name="y",
+        )
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            (1 * y) + np.ones((4, 5))
+
+    # -- matmul -----------------------------------------------------------
+
+    @pytest.mark.v1
+    @pytest.mark.parametrize("make", UNLABELED_1D)
+    def test_v1_matmul_pairs_by_size(self, xy: Variable, make: Any) -> None:
+        # matmul contracts the paired dim: length-4 array pairs with "b"
+        result = (1 * xy) @ make(range(4))
+        assert set(result.coord_dims) == {"a"}
+
+    # -- construction (add_variables bounds) ------------------------------
+
+    @pytest.mark.v1
+    def test_v1_add_variables_bound_pairs_by_size(self) -> None:
+        # The rule is the same for construction inputs: a bare-numpy bound
+        # pairs with the matching dim by size, not positionally (where the
+        # length-5 array would hit the leading dim "a" and conflict).
+        m = Model()
+        x = m.add_variables(
+            coords=[pd.RangeIndex(4, name="a"), pd.RangeIndex(5, name="time")],
+            lower=np.arange(5.0),
+            name="x",
+        )
+        assert dict(x.lower.sizes) == {"a": 4, "time": 5}
+        assert (x.lower.isel(a=0).values == np.arange(5.0)).all()
+
+    @pytest.mark.v1
+    def test_v1_add_variables_ambiguous_bound_raises(self) -> None:
+        m = Model()
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            m.add_variables(
+                coords=[pd.RangeIndex(4, name="p"), pd.RangeIndex(4, name="q")],
+                lower=np.arange(4.0),
+                name="x",
+            )
+
+    @pytest.mark.legacy
+    def test_legacy_no_divergence_is_silent(
+        self, xy: Variable, unsilenced: None
+    ) -> None:
+        # length-3 array: legacy pairs positionally with the leading dim "a";
+        # v1 would pair it with "a" too (only "a" is size 3), so there is no
+        # divergence and no warning. Pins the silent case.
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", LinopySemanticsWarning)
+            result = (1 * xy) + np.arange(3.0)
+        assert result.const.sizes == {"a": 3, "b": 4}
+        assert (result.const.isel(b=0).values == np.arange(3.0)).all()
+
+    @pytest.mark.legacy
+    def test_legacy_warns_when_v1_would_differ(
+        self, xy: Variable, unsilenced: None
+    ) -> None:
+        # length-4 array: legacy pairs positionally with "a" (size 3) → error,
+        # but the warning fires first explaining the v1 divergence.
+        # legacy positional pairing assigns the len-4 array to "a" (size 3),
+        # which then conflicts — assert it's that shape error, nothing incidental.
+        with pytest.warns(LinopySemanticsWarning, match=r"pairs by size instead"):
+            with contextlib.suppress(ValueError):
+                (1 * xy) + np.arange(4.0)
+
+    @pytest.mark.legacy
+    def test_legacy_ambiguous_pairs_positionally_with_warning(
+        self, square: Variable, unsilenced: None
+    ) -> None:
+        # The square (p:4, q:4) case where v1 *raises* — legacy must instead
+        # pair positionally with the leading dim and warn, never raise. This is
+        # the biggest legacy/v1 divergence and the strongest no-regression guard.
+        with pytest.warns(LinopySemanticsWarning, match=r"this raises"):
+            result = (1 * square) + np.arange(4.0)
+        assert result.const.sizes == {"p": 4, "q": 4}
+        # paired with the leading dim "p": the array varies along p, broadcast over q
+        assert (result.const.isel(q=0).values == np.arange(4.0)).all()
+
+    @pytest.mark.legacy
+    def test_legacy_add_variables_bound_positional(self, unsilenced: None) -> None:
+        # Regression guard for the unification: legacy add_variables still
+        # assigns an unlabeled bound positionally (here unambiguous — only "a"
+        # is size 5 — so it is also silent), producing the pre-#736 result.
+        m = Model()
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", LinopySemanticsWarning)
+            x = m.add_variables(
+                coords=[pd.RangeIndex(5, name="a"), pd.RangeIndex(3, name="b")],
+                lower=np.arange(5.0),
+                name="x",
+            )
+        assert dict(x.lower.sizes) == {"a": 5, "b": 3}
+        assert (x.lower.isel(b=0).values == np.arange(5.0)).all()
+
+
 # =====================================================================
 # §5 — User-supplied NaN raises (covers #713 and PyPSA #1683)
 # =====================================================================