diff --git a/examples/coordinate-alignment.ipynb b/examples/coordinate-alignment.ipynb
index 1547bd9d..643ab18e 100644
--- a/examples/coordinate-alignment.ipynb
+++ b/examples/coordinate-alignment.ipynb
@@ -6,13 +6,24 @@
    "source": [
     "# Coordinate Alignment\n",
     "\n",
-    "Since linopy builds on xarray, coordinate alignment matters when combining variables or expressions that live on different coordinates. By default, linopy aligns operands automatically and fills missing entries with sensible defaults. This guide shows how alignment works and how to control it with the ``join`` parameter."
+    "Since linopy builds on xarray, coordinate alignment matters when combining variables or expressions that live on different coordinates. By default, linopy uses **inner** alignment (intersection of coordinates), matching xarray's own default `arithmetic_join`. This guide shows how alignment works and how to control it with the ``join`` parameter."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.317395Z",
+     "start_time": "2026-02-19T21:28:30.389972Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:20.541634Z",
+     "iopub.status.busy": "2026-02-20T06:47:20.541299Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.421847Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.421626Z"
+    }
+   },
    "outputs": [],
    "source": [
     "import numpy as np\n",
@@ -28,13 +39,24 @@
    "source": [
     "## Default Alignment Behavior\n",
     "\n",
-    "When two operands share a dimension but have different coordinates, linopy keeps the **larger** (superset) coordinate range and fills missing positions with zeros (for addition) or zero coefficients (for multiplication)."
+    "All arithmetic operations (``+``, ``-``, ``*``, ``/``) use **\"inner\"** alignment by default — the intersection of coordinates. Only positions that exist in **both** operands appear in the result. This is safe and explicit: you always know that every entry in the result is fully defined, and nothing is silently filled in."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.379980Z",
+     "start_time": "2026-02-19T21:28:31.320575Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.423139Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.423010Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.466187Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.466002Z"
+    }
+   },
    "outputs": [],
    "source": [
     "m = linopy.Model()\n",
@@ -50,13 +72,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Adding ``x`` (5 time steps) and ``y`` (3 time steps) gives an expression over all 5 time steps. Where ``y`` has no entry (time 3, 4), the coefficient is zero — i.e. ``y`` simply drops out of the sum at those positions."
+    "Adding ``x`` (time 0–4) and ``y`` (time 0–2) gives an expression over only the **shared** time steps (0, 1, 2). Time steps 3 and 4, which only exist in ``x``, are excluded from the result."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.395989Z",
+     "start_time": "2026-02-19T21:28:31.384198Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.467497Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.467384Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.475487Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.475316Z"
+    }
+   },
    "outputs": [],
    "source": [
     "x + y"
@@ -66,13 +99,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The same applies when multiplying by a constant that covers only a subset of coordinates. Missing positions get a coefficient of zero:"
+    "The same applies when multiplying by a constant that covers only a subset of coordinates. Only the shared positions (time 0, 1, 2) appear in the result:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.411559Z",
+     "start_time": "2026-02-19T21:28:31.403124Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.476438Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.476364Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.481758Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.481585Z"
+    }
+   },
    "outputs": [],
    "source": [
     "factor = xr.DataArray([2, 3, 4], dims=[\"time\"], coords={\"time\": [0, 1, 2]})\n",
@@ -83,13 +127,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Adding a constant subset also fills missing coordinates with zero:"
+    "Adding a constant subset also restricts the result to shared coordinates:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.429941Z",
+     "start_time": "2026-02-19T21:28:31.422132Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.482751Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.482690Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.487389Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.487218Z"
+    }
+   },
    "outputs": [],
    "source": [
     "x + factor"
@@ -101,13 +156,24 @@
    "source": [
     "### Constraints with Subset RHS\n",
     "\n",
-    "For constraints, missing right-hand-side values are filled with ``NaN``, which tells linopy to **skip** the constraint at those positions:"
+    "For constraints with a DataArray right-hand side, the default alignment is ``\"left\"`` — the left operand (the expression) defines where constraints exist. Missing RHS values are filled with ``NaN``, which tells linopy to **skip** the constraint at those positions:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.441169Z",
+     "start_time": "2026-02-19T21:28:31.435220Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.488345Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.488289Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.492103Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.491931Z"
+    }
+   },
    "outputs": [],
    "source": [
     "rhs = xr.DataArray([10, 20, 30], dims=[\"time\"], coords={\"time\": [0, 1, 2]})\n",
@@ -125,53 +191,125 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "### Same-Shape Operands: Positional Alignment\n\nWhen two operands have the **same shape** on a shared dimension, linopy uses **positional alignment** by default — coordinate labels are ignored and the left operand's labels are kept. This is a performance optimization but can be surprising:"
+   "source": [
+    "### Disjoint Coordinates Produce Empty Results\n",
+    "\n",
+    "When two operands share **no** coordinate labels, the default ``\"inner\"`` alignment produces an **empty** result — the intersection of disjoint sets is empty. This makes mistakes visible immediately rather than silently filling values:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.454734Z",
+     "start_time": "2026-02-19T21:28:31.444147Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.493067Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.493009Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.499639Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.499472Z"
+    }
+   },
    "outputs": [],
    "source": [
-    "offset_const = xr.DataArray(\n",
-    "    [10, 20, 30, 40, 50], dims=[\"time\"], coords={\"time\": [5, 6, 7, 8, 9]}\n",
-    ")\n",
-    "x + offset_const"
+    "z = m.add_variables(lower=0, coords=[pd.RangeIndex(5, 10, name=\"time\")], name=\"z\")\n",
+    "x + z"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "Even though ``offset_const`` has coordinates ``[5, 6, 7, 8, 9]`` and ``x`` has ``[0, 1, 2, 3, 4]``, the result uses ``x``'s labels. The values are aligned by **position**, not by label. The same applies when adding two variables or expressions of identical shape:"
+   "source": [
+    "``x`` (time 0–4) and ``z`` (time 5–9) share no coordinate labels, so the inner join produces an **empty** expression. This is by design — if you see an empty result, it tells you the operands have no overlap, which is likely a bug or requires explicit union semantics.\n",
+    "\n",
+    "To get the union of disjoint coordinates, opt in with ``join=\"outer\"``:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.468495Z",
+     "start_time": "2026-02-19T21:28:31.461161Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.500535Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.500473Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.507040Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.506880Z"
+    }
+   },
    "outputs": [],
    "source": [
-    "z = m.add_variables(lower=0, coords=[pd.RangeIndex(5, 10, name=\"time\")], name=\"z\")\n",
-    "x + z"
+    "x.add(z, join=\"outer\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "``x`` (time 0–4) and ``z`` (time 5–9) share no coordinate labels, yet the result has 5 entries under ``x``'s coordinates — because they have the same shape, positions are matched directly.\n\nTo force **label-based** alignment, pass an explicit ``join``:"
+   "source": [
+    "With ``join=\"outer\"``, the result spans all 10 time steps. Each variable appears only at its own positions; missing entries are filled with zero coefficients.\n",
+    "\n",
+    "The same works for constant operands with different labels:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:35:16.953353Z",
+     "start_time": "2026-02-19T21:35:16.927666Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.507911Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.507854Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.512952Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.512792Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "offset_const = xr.DataArray(\n",
+    "    [10, 20, 30, 40, 50], dims=[\"time\"], coords={\"time\": [5, 6, 7, 8, 9]}\n",
+    ")\n",
+    "x.add(offset_const, join=\"outer\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
    "metadata": {},
+   "source": [
+    "The result spans 10 time steps. At time 0–4 only ``x`` contributes; at time 5–9 only the constant applies. Missing entries are filled with appropriate defaults (zero for addition).\n",
+    "\n",
+    "To force **positional** alignment (ignoring labels), use ``join=\"override\"``:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.513886Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.513834Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.519693Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.519438Z"
+    }
+   },
    "outputs": [],
    "source": [
-    "x.add(z, join=\"outer\")"
+    "x.add(z, join=\"override\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "With ``join=\"outer\"``, the result spans all 10 time steps (union of 0–4 and 5–9), filling missing positions with zeros. This is the correct label-based alignment. The same-shape positional shortcut is equivalent to ``join=\"override\"`` — see below."
+   "source": [
+    "With ``join=\"override\"``, positions are matched directly and the left operand's labels are used. This can be useful when coordinate labels are irrelevant and you just want to pair entries by position."
+   ]
   },
   {
    "cell_type": "markdown",
@@ -181,18 +319,31 @@
     "\n",
     "For explicit control over alignment, use the ``.add()``, ``.sub()``, ``.mul()``, and ``.div()`` methods with a ``join`` parameter. The supported values follow xarray conventions:\n",
     "\n",
-    "- ``\"inner\"`` — intersection of coordinates\n",
-    "- ``\"outer\"`` — union of coordinates (with fill)\n",
+    "- ``\"inner\"`` (default) — intersection of coordinates\n",
+    "- ``\"outer\"`` — union of coordinates (with fill values)\n",
     "- ``\"left\"`` — keep left operand's coordinates\n",
     "- ``\"right\"`` — keep right operand's coordinates\n",
     "- ``\"override\"`` — positional alignment, ignore coordinate labels\n",
-    "- ``\"exact\"`` — coordinates must match exactly (raises on mismatch)"
+    "- ``\"exact\"`` — coordinates must match exactly (raises on mismatch)\n",
+    "\n",
+    "Note: operators (``+``, ``-``, ``*``, ``/``) always use the default ``\"inner\"`` join. Use the named methods to specify a different join."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.495689Z",
+     "start_time": "2026-02-19T21:28:31.490736Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.520655Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.520600Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.524087Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.523920Z"
+    }
+   },
    "outputs": [],
    "source": [
     "m2 = linopy.Model()\n",
@@ -208,13 +359,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "**Inner join** — only shared coordinates (i=1, 2):"
+    "**Inner join** (the default) — only shared coordinates (i=1, 2). This is what ``a + b`` would produce:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.511734Z",
+     "start_time": "2026-02-19T21:28:31.502665Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.524955Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.524901Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.531133Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.530941Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.add(b, join=\"inner\")"
@@ -230,7 +392,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.573137Z",
+     "start_time": "2026-02-19T21:28:31.565103Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.532103Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.532044Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.538610Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.538441Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.add(b, join=\"outer\")"
@@ -246,7 +419,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.589234Z",
+     "start_time": "2026-02-19T21:28:31.580295Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.539516Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.539459Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.545958Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.545770Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.add(b, join=\"left\")"
@@ -262,7 +446,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.605303Z",
+     "start_time": "2026-02-19T21:28:31.597026Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.546895Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.546835Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.553932Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.553719Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.add(b, join=\"right\")"
@@ -271,12 +466,25 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "**Override** — positional alignment, ignore coordinate labels. The result uses the left operand's coordinates. Here ``a`` has i=[0, 1, 2] and ``b`` has i=[1, 2, 3], so positions are matched as 0↔1, 1↔2, 2↔3:"
+   "source": [
+    "**Override** — positional alignment, ignore coordinate labels. The result uses the left operand's coordinates. Here ``a`` has i=[0, 1, 2] and ``b`` has i=[1, 2, 3], so positions are matched as 0↔1, 1↔2, 2↔3:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.622016Z",
+     "start_time": "2026-02-19T21:28:31.613412Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.555055Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.554996Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.561786Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.561201Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.add(b, join=\"override\")"
@@ -288,13 +496,24 @@
    "source": [
     "### Multiplication with ``join``\n",
     "\n",
-    "The same ``join`` parameter works on ``.mul()`` and ``.div()``. When multiplying by a constant that covers a subset, ``join=\"inner\"`` restricts the result to shared coordinates only, while ``join=\"left\"`` fills missing values with zero:"
+    "The same ``join`` parameter works on ``.mul()`` and ``.div()``. When multiplying by a constant that covers a subset, ``join=\"inner\"`` (the default) restricts the result to shared coordinates only, while ``join=\"left\"`` keeps the left operand's coordinates and fills missing values with zero:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.634088Z",
+     "start_time": "2026-02-19T21:28:31.627066Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.562787Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.562721Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.568020Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.567846Z"
+    }
+   },
    "outputs": [],
    "source": [
     "const = xr.DataArray([2, 3, 4], dims=[\"i\"], coords={\"i\": [1, 2, 3]})\n",
@@ -305,7 +524,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.654190Z",
+     "start_time": "2026-02-19T21:28:31.648732Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.568990Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.568925Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.573133Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.572973Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.mul(const, join=\"left\")"
@@ -317,13 +547,26 @@
    "source": [
     "## Alignment in Constraints\n",
     "\n",
-    "The ``.le()``, ``.ge()``, and ``.eq()`` methods create constraints with explicit coordinate alignment. They accept the same ``join`` parameter:"
+    "For constraints with a DataArray right-hand side, the default alignment is ``\"left\"`` — the expression defines where constraints exist, and missing RHS values become ``NaN`` (meaning \"no constraint\").\n",
+    "\n",
+    "The ``.le()``, ``.ge()``, and ``.eq()`` methods accept the same ``join`` parameter for explicit control:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.675417Z",
+     "start_time": "2026-02-19T21:28:31.668026Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.574043Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.573985Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.578196Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.578028Z"
+    }
+   },
    "outputs": [],
    "source": [
     "rhs = xr.DataArray([10, 20], dims=[\"i\"], coords={\"i\": [0, 1]})\n",
@@ -335,13 +578,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "With ``join=\"inner\"``, the constraint only exists at the intersection (i=0, 1). Compare with ``join=\"left\"``:"
+    "With ``join=\"inner\"``, the constraint only exists at the intersection (i=0, 1). Compare with the default ``join=\"left\"``:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.690980Z",
+     "start_time": "2026-02-19T21:28:31.685057Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.579086Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.579030Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.582636Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.582451Z"
+    }
+   },
    "outputs": [],
    "source": [
     "a.le(rhs, join=\"left\")"
@@ -351,7 +605,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "With ``join=\"left\"``, the result covers all of ``a``'s coordinates (i=0, 1, 2). At i=2, where the RHS has no value, the RHS becomes ``NaN`` and the constraint is masked out.\n",
+    "With ``join=\"left\"`` (the default for constraints with DataArray RHS), the result covers all of ``a``'s coordinates (i=0, 1, 2). At i=2, where the RHS has no value, the RHS becomes ``NaN`` and the constraint is masked out.\n",
     "\n",
     "The same methods work on expressions:"
    ]
@@ -359,7 +613,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.703444Z",
+     "start_time": "2026-02-19T21:28:31.694875Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.583521Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.583468Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.590082Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.589917Z"
+    }
+   },
    "outputs": [],
    "source": [
     "expr = 2 * a + 1\n",
@@ -369,12 +634,27 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Practical Example\n\nConsider a generation dispatch model where solar availability follows a daily profile and a minimum demand constraint only applies during peak hours."
+   "source": [
+    "## Practical Example\n",
+    "\n",
+    "Consider a generation dispatch model where solar availability follows a daily profile and a minimum demand constraint only applies during peak hours."
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.714348Z",
+     "start_time": "2026-02-19T21:28:31.709096Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.590977Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.590921Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.594137Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.593957Z"
+    }
+   },
    "outputs": [],
    "source": [
     "m3 = linopy.Model()\n",
@@ -395,7 +675,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.734205Z",
+     "start_time": "2026-02-19T21:28:31.721553Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.595024Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.594969Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.604003Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.603809Z"
+    }
+   },
    "outputs": [],
    "source": [
     "capacity = xr.DataArray([100, 80, 50], dims=[\"tech\"], coords={\"tech\": techs})\n",
@@ -405,12 +696,25 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "For solar, we build a full 24-hour availability profile — zero at night, sine-shaped during daylight (hours 6–18). Since this covers all hours, standard alignment works directly and solar is properly constrained to zero at night:"
+   "source": [
+    "For solar, we build a full 24-hour availability profile — zero at night, sine-shaped during daylight (hours 6–18). Since this covers all hours, standard alignment works directly and solar is properly constrained to zero at night:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.753434Z",
+     "start_time": "2026-02-19T21:28:31.741594Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.604903Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.604846Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.613203Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.613030Z"
+    }
+   },
    "outputs": [],
    "source": [
     "solar_avail = np.zeros(24)\n",
@@ -424,12 +728,25 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "Now suppose a minimum demand of 120 MW must be met, but only during peak hours (8–20). The demand array covers a subset of hours, so we use ``join=\"inner\"`` to restrict the constraint to just those hours:"
+   "source": [
+    "Now suppose a minimum demand of 120 MW must be met, but only during peak hours (8–20). The demand array covers a subset of hours, so we use ``join=\"inner\"`` to restrict the constraint to just those hours:"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-02-19T21:28:31.768893Z",
+     "start_time": "2026-02-19T21:28:31.756306Z"
+    },
+    "execution": {
+     "iopub.execute_input": "2026-02-20T06:47:21.614093Z",
+     "iopub.status.busy": "2026-02-20T06:47:21.614038Z",
+     "iopub.status.idle": "2026-02-20T06:47:21.623564Z",
+     "shell.execute_reply": "2026-02-20T06:47:21.623388Z"
+    }
+   },
    "outputs": [],
    "source": [
     "peak_hours = pd.RangeIndex(8, 21, name=\"hour\")\n",
@@ -444,7 +761,9 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "The demand constraint only applies during peak hours (8–20). Outside that range, no minimum generation is required."
+   "source": [
+    "The demand constraint only applies during peak hours (8–20). Outside that range, no minimum generation is required."
+   ]
   },
   {
    "cell_type": "markdown",
@@ -452,12 +771,17 @@
    "source": [
     "## Summary\n",
     "\n",
+    "| Context | Default ``join`` | Behavior |\n",
+    "|---------|-----------------|----------|\n",
+    "| Arithmetic (``+``, ``-``, ``*``, ``/``) | ``\"inner\"`` | Intersection of coordinates; safe and explicit |\n",
+    "| Constraint with DataArray RHS | ``\"left\"`` | Expression defines positions; missing RHS becomes NaN (masked) |\n",
+    "| Constraint with expression RHS | ``\"inner\"`` | Goes through subtraction, inherits arithmetic default |\n",
+    "\n",
     "| ``join`` | Coordinates | Fill behavior |\n",
     "|----------|------------|---------------|\n",
-    "| ``None`` (default) | Auto-detect (keeps superset) | Zeros for arithmetic, NaN for constraint RHS |\n",
-    "| ``\"inner\"`` | Intersection only | No fill needed |\n",
-    "| ``\"outer\"`` | Union | Fill with operation identity (0 for add, 0 for mul) |\n",
-    "| ``\"left\"`` | Left operand's | Fill right with identity |\n",
+    "| ``\"inner\"`` (arithmetic default) | Intersection only | No fill needed |\n",
+    "| ``\"outer\"`` | Union | Fill with operation identity (0 for add, 0 for mul, 1 for div) |\n",
+    "| ``\"left\"`` (constraint default) | Left operand's | Fill right with identity / NaN for constraint |\n",
     "| ``\"right\"`` | Right operand's | Fill left with identity |\n",
     "| ``\"override\"`` | Left operand's (positional) | Positional alignment, ignore labels |\n",
     "| ``\"exact\"`` | Must match exactly | Raises error if different |"
@@ -480,7 +804,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,
diff --git a/linopy/expressions.py b/linopy/expressions.py
index e1fbe1a9..24714f00 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -48,7 +48,6 @@
     LocIndexer,
     as_dataarray,
     assign_multiindex_safe,
-    check_common_keys_values,
     check_has_nulls,
     check_has_nulls_polars,
     fill_missing_coords,
@@ -537,9 +536,9 @@ def _align_constant(
         other : DataArray
             The constant to align.
         fill_value : float, default: 0
-            Fill value for missing coordinates.
+            Fill value for missing coordinates in the other operand.
         join : str, optional
-            Alignment method. If None, uses size-aware default behavior.
+            Alignment method. If None, defaults to "outer".
 
         Returns
         -------
@@ -551,19 +550,20 @@ def _align_constant(
             Whether the expression's data needs reindexing.
         """
         if join is None:
-            if other.sizes == self.const.sizes:
-                return self.const, other.assign_coords(coords=self.coords), False
+            join = "inner"
+
+        if join == "override":
+            return self.const, other.assign_coords(coords=self.coords), False
+        elif join == "left":
             return (
                 self.const,
                 other.reindex_like(self.const, fill_value=fill_value),
                 False,
             )
-        elif join == "override":
-            return self.const, other.assign_coords(coords=self.coords), False
         else:
-            self_const, aligned = xr.align(
-                self.const, other, join=join, fill_value=fill_value
-            )
+            self_const, aligned = xr.align(self.const, other, join=join)
+            self_const = self_const.fillna(0)  # expression const: always 0 when missing
+            aligned = aligned.fillna(fill_value)  # operand fill: operation-dependent
             return self_const, aligned, True
 
     def _add_constant(
@@ -597,9 +597,12 @@ def _apply_constant_op(
         )
         if needs_data_reindex:
             data = self.data.reindex_like(self_const, fill_value=self._fill_value)
+            # Positions where expression was missing get coeffs=NaN from fill.
+            # Replace with 0 so op(0, factor) = 0 instead of NaN.
+            coeffs = data.coeffs.fillna(0)
             return self.__class__(
                 assign_multiindex_safe(
-                    data, coeffs=op(data.coeffs, factor), const=op(self_const, factor)
+                    data, coeffs=op(coeffs, factor), const=op(self_const, factor)
                 ),
                 self.model,
             )
@@ -608,7 +611,7 @@ def _apply_constant_op(
     def _multiply_by_constant(
         self: GenericExpression, other: ConstantLike, join: str | None = None
     ) -> GenericExpression:
-        return self._apply_constant_op(other, operator.mul, fill_value=0, join=join)
+        return self._apply_constant_op(other, operator.mul, fill_value=1, join=join)
 
     def _divide_by_constant(
         self: GenericExpression, other: ConstantLike, join: str | None = None
@@ -1082,8 +1085,31 @@ def to_constraint(
                     f"RHS DataArray has dimensions {extra_dims} not present "
                     f"in the expression. Cannot create constraint."
                 )
-            rhs = rhs.reindex_like(self.const, fill_value=np.nan)
-        elif isinstance(rhs, np.ndarray | pd.Series | pd.DataFrame) and rhs.ndim > len(
+            # Align rhs to expression coordinates. NaN = "no constraint".
+            if join is None or join == "left":
+                aligned_rhs = rhs.reindex_like(self.const, fill_value=np.nan)
+                expr_data = self.data
+                expr_const = self.const
+            elif join == "override":
+                aligned_rhs = rhs.assign_coords(coords=self.const.coords)
+                expr_data = self.data
+                expr_const = self.const
+            else:
+                expr_const_aligned, aligned_rhs = xr.align(
+                    self.const, rhs, join=join, fill_value=np.nan
+                )
+                expr_data = self.data.reindex_like(
+                    expr_const_aligned, fill_value=self._fill_value
+                )
+                expr_const = expr_const_aligned.fillna(0)
+            # Compute constraint directly: rhs = aligned_rhs - expr_const
+            constraint_rhs = aligned_rhs - expr_const
+            data = assign_multiindex_safe(
+                expr_data[["coeffs", "vars"]], sign=sign, rhs=constraint_rhs
+            )
+            return constraints.Constraint(data, model=self.model)
+
+        if isinstance(rhs, np.ndarray | pd.Series | pd.DataFrame) and rhs.ndim > len(
             self.coord_dims
         ):
             raise ValueError(
@@ -2320,16 +2346,6 @@ def merge(
 
     model = exprs[0].model
 
-    if join is not None:
-        override = join == "override"
-    elif cls in linopy_types and dim in HELPER_DIMS:
-        coord_dims = [
-            {k: v for k, v in e.sizes.items() if k not in HELPER_DIMS} for e in exprs
-        ]
-        override = check_common_keys_values(coord_dims)  # type: ignore
-    else:
-        override = False
-
     data = [e.data if isinstance(e, linopy_types) else e for e in exprs]
     data = [fill_missing_coords(ds, fill_helper_dims=True) for ds in data]
 
@@ -2345,10 +2361,8 @@ def merge(
 
         if join is not None:
             kwargs["join"] = join
-        elif override:
-            kwargs["join"] = "override"
         else:
-            kwargs.setdefault("join", "outer")
+            kwargs["join"] = "inner"
 
     if dim == TERM_DIM:
         ds = xr.concat([d[["coeffs", "vars"]] for d in data], dim, **kwargs)
@@ -2361,6 +2375,20 @@ def merge(
         const = xr.concat([d["const"] for d in data], dim, **kwargs).prod(FACTOR_DIM)
         ds = assign_multiindex_safe(ds, coeffs=coeffs, const=const)
     else:
+        # When concatenating along a coordinate dimension, pre-pad helper
+        # dimensions (_term, _factor) to the same size so the join setting
+        # only affects coordinate dimensions (not the term/factor axes).
+        fill = kwargs.get("fill_value", FILL_VALUE)
+        for helper_dim in HELPER_DIMS:
+            sizes = [d.sizes.get(helper_dim, 0) for d in data]
+            max_size = max(sizes) if sizes else 0
+            if max_size > 0 and min(sizes) < max_size:
+                data = [
+                    d.reindex({helper_dim: range(max_size)}, fill_value=fill)
+                    if d.sizes.get(helper_dim, 0) < max_size
+                    else d
+                    for d in data
+                ]
         ds = xr.concat(data, dim, **kwargs)
 
     for d in set(HELPER_DIMS) & set(ds.coords):
diff --git a/linopy/variables.py b/linopy/variables.py
index 0eea6634..cf137f8c 100644
--- a/linopy/variables.py
+++ b/linopy/variables.py
@@ -401,7 +401,10 @@ def __mul__(self, other: SideLike) -> ExpressionLike:
             if isinstance(other, Variable | ScalarVariable):
                 return self.to_linexpr() * other
 
-            return self.to_linexpr(other)
+            if isinstance(other, expressions.LinearExpression):
+                return self.to_linexpr() * other
+
+            return self.to_linexpr()._multiply_by_constant(other)
         except TypeError:
             return NotImplemented
 
diff --git a/test/test_linear_expression.py b/test/test_linear_expression.py
index 2af1a8ea..c1fa5f9b 100644
--- a/test/test_linear_expression.py
+++ b/test/test_linear_expression.py
@@ -443,10 +443,14 @@ def test_linear_expression_sum(
 
     assert_linequal(expr.sum(["dim_0", TERM_DIM]), expr.sum("dim_0"))
 
-    # test special case otherride coords
+    # With "inner" default, disjoint coord slices produce empty intersection
     expr = v.loc[:9] + v.loc[10:]
+    assert len(expr.coords["dim_2"]) == 0
+
+    # Use join="outer" to get the full union
+    expr = v.loc[:9].add(v.loc[10:], join="outer")
     assert expr.nterm == 2
-    assert len(expr.coords["dim_2"]) == 10
+    assert len(expr.coords["dim_2"]) == 20
 
 
 def test_linear_expression_sum_with_const(
@@ -467,8 +471,8 @@ def test_linear_expression_sum_with_const(
 
     assert_linequal(expr.sum(["dim_0", TERM_DIM]), expr.sum("dim_0"))
 
-    # test special case otherride coords
-    expr = v.loc[:9] + v.loc[10:]
+    # With "outer" default, disjoint coord slices produce the full union
+    expr = v.loc[:9].add(v.loc[10:], "override")
     assert expr.nterm == 2
     assert len(expr.coords["dim_2"]) == 10
 
@@ -590,16 +594,14 @@ def superset(self) -> xr.DataArray:
 
     @pytest.fixture
     def expected_fill(self) -> np.ndarray:
-        arr = np.zeros(20)
-        arr[1] = 10.0
-        arr[3] = 30.0
-        return arr
+        """Expected values at intersection coords [1, 3]."""
+        return np.array([10.0, 30.0])
 
     def test_var_mul_subset(
         self, v: Variable, subset: xr.DataArray, expected_fill: np.ndarray
     ) -> None:
         result = v * subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         np.testing.assert_array_equal(result.coeffs.squeeze().values, expected_fill)
 
@@ -608,7 +610,7 @@ def test_expr_mul_subset(
     ) -> None:
         expr = 1 * v
         result = expr * subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         np.testing.assert_array_equal(result.coeffs.squeeze().values, expected_fill)
 
@@ -631,7 +633,7 @@ def test_var_add_subset(
         self, v: Variable, subset: xr.DataArray, expected_fill: np.ndarray
     ) -> None:
         result = v + subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, expected_fill)
 
@@ -639,7 +641,7 @@ def test_var_sub_subset(
         self, v: Variable, subset: xr.DataArray, expected_fill: np.ndarray
     ) -> None:
         result = v - subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, -expected_fill)
 
@@ -651,7 +653,7 @@ def test_expr_add_subset(
     ) -> None:
         expr = v + 5
         result = expr + subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, expected_fill + 5)
 
@@ -660,7 +662,7 @@ def test_expr_sub_subset(
     ) -> None:
         expr = v + 5
         result = expr - subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, 5 - expected_fill)
 
@@ -670,10 +672,10 @@ def test_subset_sub_expr(self, v: Variable, subset: xr.DataArray) -> None:
 
     def test_var_div_subset(self, v: Variable, subset: xr.DataArray) -> None:
         result = v / subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         assert result.coeffs.squeeze().sel(dim_2=1).item() == pytest.approx(0.1)
-        assert result.coeffs.squeeze().sel(dim_2=0).item() == pytest.approx(1.0)
+        assert result.coeffs.squeeze().sel(dim_2=3).item() == pytest.approx(1 / 30)
 
     def test_var_le_subset(self, v: Variable, subset: xr.DataArray) -> None:
         con = v <= subset
@@ -712,11 +714,11 @@ def test_add_commutativity_full_coords(self, v: Variable) -> None:
         )
         assert_linequal(v + full, full + v)
 
-    def test_superset_addition_pins_to_lhs(
+    def test_superset_addition_keeps_intersection(
         self, v: Variable, superset: xr.DataArray
     ) -> None:
         result = v + superset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == v.sizes["dim_2"]  # inner = v's coords
         assert not np.isnan(result.const.values).any()
 
     def test_superset_add_var(self, v: Variable, superset: xr.DataArray) -> None:
@@ -740,34 +742,34 @@ def test_superset_comparison_var(
         assert not np.isnan(con.lhs.coeffs.values).any()
         assert not np.isnan(con.rhs.values).any()
 
-    def test_disjoint_addition_pins_to_lhs(self, v: Variable) -> None:
+    def test_disjoint_addition_gives_empty(self, v: Variable) -> None:
         disjoint = xr.DataArray(
             [100.0, 200.0], dims=["dim_2"], coords={"dim_2": [50, 60]}
         )
         result = v + disjoint
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
-        assert not np.isnan(result.const.values).any()
-        np.testing.assert_array_equal(result.const.values, np.zeros(20))
+        assert result.sizes["dim_2"] == 0  # inner: no shared coords
 
     def test_expr_div_subset(self, v: Variable, subset: xr.DataArray) -> None:
         expr = 1 * v
         result = expr / subset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         assert result.coeffs.squeeze().sel(dim_2=1).item() == pytest.approx(0.1)
-        assert result.coeffs.squeeze().sel(dim_2=0).item() == pytest.approx(1.0)
+        assert result.coeffs.squeeze().sel(dim_2=3).item() == pytest.approx(1 / 30)
 
     def test_subset_add_var_coefficients(
         self, v: Variable, subset: xr.DataArray
     ) -> None:
         result = subset + v
-        np.testing.assert_array_equal(result.coeffs.squeeze().values, np.ones(20))
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
+        np.testing.assert_array_equal(result.coeffs.squeeze().values, np.ones(2))
 
     def test_subset_sub_var_coefficients(
         self, v: Variable, subset: xr.DataArray
     ) -> None:
         result = subset - v
-        np.testing.assert_array_equal(result.coeffs.squeeze().values, -np.ones(20))
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
+        np.testing.assert_array_equal(result.coeffs.squeeze().values, -np.ones(2))
 
     @pytest.mark.parametrize("sign", ["<=", ">=", "=="])
     def test_subset_comparison_var(
@@ -783,21 +785,21 @@ def test_subset_comparison_var(
         assert np.isnan(con.rhs.sel(dim_2=0).item())
         assert con.rhs.sel(dim_2=1).item() == pytest.approx(10.0)
 
-    def test_superset_mul_pins_to_lhs(
+    def test_superset_mul_keeps_intersection(
         self, v: Variable, superset: xr.DataArray
     ) -> None:
         result = v * superset
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == v.sizes["dim_2"]  # inner = v's coords
         assert not np.isnan(result.coeffs.values).any()
 
-    def test_superset_div_pins_to_lhs(self, v: Variable) -> None:
+    def test_superset_div_keeps_intersection(self, v: Variable) -> None:
         superset_nonzero = xr.DataArray(
             np.arange(1, 26, dtype=float),
             dims=["dim_2"],
             coords={"dim_2": range(25)},
         )
         result = v / superset_nonzero
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == v.sizes["dim_2"]  # inner = v's coords
         assert not np.isnan(result.coeffs.values).any()
 
     def test_quadexpr_add_subset(
@@ -806,7 +808,7 @@ def test_quadexpr_add_subset(
         qexpr = v * v
         result = qexpr + subset
         assert isinstance(result, QuadraticExpression)
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, expected_fill)
 
@@ -816,7 +818,7 @@ def test_quadexpr_sub_subset(
         qexpr = v * v
         result = qexpr - subset
         assert isinstance(result, QuadraticExpression)
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.const.values).any()
         np.testing.assert_array_equal(result.const.values, -expected_fill)
 
@@ -826,7 +828,7 @@ def test_quadexpr_mul_subset(
         qexpr = v * v
         result = qexpr * subset
         assert isinstance(result, QuadraticExpression)
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         np.testing.assert_array_equal(result.coeffs.squeeze().values, expected_fill)
 
@@ -836,7 +838,7 @@ def test_subset_mul_quadexpr(
         qexpr = v * v
         result = subset * qexpr
         assert isinstance(result, QuadraticExpression)
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
+        assert result.sizes["dim_2"] == 2  # inner: intersection [1, 3]
         assert not np.isnan(result.coeffs.values).any()
         np.testing.assert_array_equal(result.coeffs.squeeze().values, expected_fill)
 
@@ -855,13 +857,11 @@ def test_multidim_subset_mul(self, m: Model) -> None:
             coords={"a": [1, 3], "b": [0, 4]},
         )
         result = w * subset_2d
-        assert result.sizes["a"] == 4
-        assert result.sizes["b"] == 5
+        assert result.sizes["a"] == 2  # inner: intersection [1, 3]
+        assert result.sizes["b"] == 2  # inner: intersection [0, 4]
         assert not np.isnan(result.coeffs.values).any()
         assert result.coeffs.squeeze().sel(a=1, b=0).item() == pytest.approx(2.0)
         assert result.coeffs.squeeze().sel(a=3, b=4).item() == pytest.approx(5.0)
-        assert result.coeffs.squeeze().sel(a=0, b=0).item() == pytest.approx(0.0)
-        assert result.coeffs.squeeze().sel(a=1, b=2).item() == pytest.approx(0.0)
 
     def test_multidim_subset_add(self, m: Model) -> None:
         coords_a = pd.RangeIndex(4, name="a")
@@ -874,12 +874,11 @@ def test_multidim_subset_add(self, m: Model) -> None:
             coords={"a": [1, 3], "b": [0, 4]},
         )
         result = w + subset_2d
-        assert result.sizes["a"] == 4
-        assert result.sizes["b"] == 5
+        assert result.sizes["a"] == 2  # inner: intersection [1, 3]
+        assert result.sizes["b"] == 2  # inner: intersection [0, 4]
         assert not np.isnan(result.const.values).any()
         assert result.const.sel(a=1, b=0).item() == pytest.approx(2.0)
         assert result.const.sel(a=3, b=4).item() == pytest.approx(5.0)
-        assert result.const.sel(a=0, b=0).item() == pytest.approx(0.0)
 
     def test_constraint_rhs_extra_dims_raises(self, v: Variable) -> None:
         rhs = xr.DataArray(
@@ -893,23 +892,19 @@ def test_da_truediv_var_raises(self, v: Variable) -> None:
         with pytest.raises(TypeError):
             da / v  # type: ignore[operator]
 
-    def test_disjoint_mul_produces_zeros(self, v: Variable) -> None:
+    def test_disjoint_mul_gives_empty(self, v: Variable) -> None:
         disjoint = xr.DataArray(
             [10.0, 20.0], dims=["dim_2"], coords={"dim_2": [50, 60]}
         )
         result = v * disjoint
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
-        assert not np.isnan(result.coeffs.values).any()
-        np.testing.assert_array_equal(result.coeffs.squeeze().values, np.zeros(20))
+        assert result.sizes["dim_2"] == 0  # inner: no shared coords
 
-    def test_disjoint_div_preserves_coeffs(self, v: Variable) -> None:
+    def test_disjoint_div_gives_empty(self, v: Variable) -> None:
         disjoint = xr.DataArray(
             [10.0, 20.0], dims=["dim_2"], coords={"dim_2": [50, 60]}
         )
         result = v / disjoint
-        assert result.sizes["dim_2"] == v.sizes["dim_2"]
-        assert not np.isnan(result.coeffs.values).any()
-        np.testing.assert_array_equal(result.coeffs.squeeze().values, np.ones(20))
+        assert result.sizes["dim_2"] == 0  # inner: no shared coords
 
     def test_da_eq_da_still_works(self) -> None:
         da1 = xr.DataArray([1, 2, 3])
@@ -1840,7 +1835,7 @@ def test_mul_constant_join_outer(self, a: Variable) -> None:
         const = xr.DataArray([2, 3, 4], dims=["i"], coords={"i": [1, 2, 3]})
         result = a.to_linexpr().mul(const, join="outer")
         assert list(result.data.indexes["i"]) == [0, 1, 2, 3]
-        assert result.coeffs.sel(i=0).item() == 0
+        assert result.coeffs.sel(i=0).item() == 1  # fill=1 (identity for mul)
         assert result.coeffs.sel(i=1).item() == 2
         assert result.coeffs.sel(i=2).item() == 3
 
@@ -1956,12 +1951,13 @@ def test_mul_constant_outer_fill_values(self, a: Variable) -> None:
         other = xr.DataArray([2, 3], dims=["i"], coords={"i": [1, 3]})
         result = expr.mul(other, join="outer")
         assert set(result.coords["i"].values) == {0, 1, 2, 3}
-        assert result.const.sel(i=0).item() == 0
-        assert result.const.sel(i=1).item() == 10
-        assert result.const.sel(i=2).item() == 0
-        assert result.const.sel(i=3).item() == 0
+        # fill=1: missing factor positions keep expression unchanged
+        assert result.const.sel(i=0).item() == 5  # 5 * 1
+        assert result.const.sel(i=1).item() == 10  # 5 * 2
+        assert result.const.sel(i=2).item() == 5  # 5 * 1
+        assert result.const.sel(i=3).item() == 0  # 0 * 3 (no expression here)
         assert result.coeffs.squeeze().sel(i=1).item() == 2
-        assert result.coeffs.squeeze().sel(i=0).item() == 0
+        assert result.coeffs.squeeze().sel(i=0).item() == 1  # 1 * 1 (fill=1)
 
     def test_div_constant_override_positional(self, a: Variable) -> None:
         expr = 1 * a + 10
@@ -2025,10 +2021,10 @@ def test_add_scalar_with_explicit_join(self, a: Variable) -> None:
         assert list(result.coords["i"].values) == [0, 1, 2]
 
     def test_quadratic_add_constant_join_inner(self, a: Variable, b: Variable) -> None:
-        quad = a.to_linexpr() * b.to_linexpr()
+        quad = a.to_linexpr() * b.to_linexpr()  # inner default → i=[1, 2]
         const = xr.DataArray([10, 20, 30], dims=["i"], coords={"i": [1, 2, 3]})
         result = quad.add(const, join="inner")
-        assert list(result.data.indexes["i"]) == [1, 2, 3]
+        assert list(result.data.indexes["i"]) == [1, 2]
 
     def test_quadratic_add_expr_join_inner(self, a: Variable) -> None:
         quad = a.to_linexpr() * a.to_linexpr()
@@ -2037,10 +2033,10 @@ def test_quadratic_add_expr_join_inner(self, a: Variable) -> None:
         assert list(result.data.indexes["i"]) == [0, 1]
 
     def test_quadratic_mul_constant_join_inner(self, a: Variable, b: Variable) -> None:
-        quad = a.to_linexpr() * b.to_linexpr()
+        quad = a.to_linexpr() * b.to_linexpr()  # inner default → i=[1, 2]
         const = xr.DataArray([2, 3, 4], dims=["i"], coords={"i": [1, 2, 3]})
         result = quad.mul(const, join="inner")
-        assert list(result.data.indexes["i"]) == [1, 2, 3]
+        assert list(result.data.indexes["i"]) == [1, 2]
 
     def test_merge_join_left(self, a: Variable, b: Variable) -> None:
         result: LinearExpression = merge([a.to_linexpr(), b.to_linexpr()], join="left")
@@ -2049,3 +2045,79 @@ def test_merge_join_left(self, a: Variable, b: Variable) -> None:
     def test_merge_join_right(self, a: Variable, b: Variable) -> None:
         result: LinearExpression = merge([a.to_linexpr(), b.to_linexpr()], join="right")
         assert list(result.data.indexes["i"]) == [1, 2, 3]
+
+    def test_constraint_rhs_join_inner(self, a: Variable) -> None:
+        rhs = xr.DataArray([10, 20], dims=["i"], coords={"i": [1, 2]})
+        con = a.to_linexpr().le(rhs, join="inner")
+        assert list(con.data.indexes["i"]) == [1, 2]
+        assert con.rhs.sel(i=1).item() == 10.0
+        assert con.rhs.sel(i=2).item() == 20.0
+
+    def test_constraint_rhs_join_outer(self, a: Variable) -> None:
+        rhs = xr.DataArray([10, 20], dims=["i"], coords={"i": [1, 2]})
+        con = a.to_linexpr().le(rhs, join="outer")
+        assert list(con.data.indexes["i"]) == [0, 1, 2]
+        assert np.isnan(con.rhs.sel(i=0).item())
+        assert con.rhs.sel(i=1).item() == 10.0
+        assert con.rhs.sel(i=2).item() == 20.0
+
+
+class TestAssociativity:
+    """Verify that addition is associative with the 'inner' default."""
+
+    @pytest.fixture
+    def m3(self) -> Model:
+        m = Model()
+        m.add_variables(coords=[pd.Index([0, 1, 2], name="t")], name="y")
+        m.add_variables(coords=[pd.Index(range(5), name="t")], name="x")
+        return m
+
+    def test_associativity_with_constant(self, m3: Model) -> None:
+        y = m3.variables["y"]
+        x = m3.variables["x"]
+        factor = xr.DataArray(
+            np.arange(6, dtype=float), dims=["t"], coords={"t": range(6)}
+        )
+        r1 = (y + x) + factor
+        r2 = (y + factor) + x
+        r3 = y + (x + factor)
+        assert_linequal(r1, r2)
+        assert_linequal(r1, r3)
+        # inner: intersection of all three is t=[0,1,2]
+        assert set(r1.coords["t"].values) == {0, 1, 2}
+
+    def test_same_shape_disjoint_add(self, m3: Model) -> None:
+        m = Model()
+        a = m.add_variables(coords=[pd.Index([0, 1, 2], name="i")], name="a")
+        b = m.add_variables(coords=[pd.Index([5, 6, 7], name="i")], name="b")
+        result = a + b
+        assert result.sizes["i"] == 0  # inner: no shared coords
+
+    def test_same_shape_disjoint_mul(self, m3: Model) -> None:
+        m = Model()
+        a = m.add_variables(coords=[pd.Index([0, 1, 2], name="i")], name="a")
+        factor = xr.DataArray([10, 20, 30], dims=["i"], coords={"i": [5, 6, 7]})
+        result = a * factor
+        assert result.sizes["i"] == 0  # inner: no shared coords
+
+    def test_same_shape_disjoint_div(self, m3: Model) -> None:
+        m = Model()
+        a = m.add_variables(coords=[pd.Index([0, 1, 2], name="i")], name="a")
+        divisor = xr.DataArray([10.0, 20.0, 30.0], dims=["i"], coords={"i": [5, 6, 7]})
+        result = a / divisor
+        assert result.sizes["i"] == 0  # inner: no shared coords
+
+    def test_override_gives_positional_matching(self, m3: Model) -> None:
+        m = Model()
+        a = m.add_variables(coords=[pd.Index([0, 1, 2], name="i")], name="a")
+        b = m.add_variables(coords=[pd.Index([5, 6, 7], name="i")], name="b")
+        result = a.add(b, join="override")
+        assert list(result.coords["i"].values) == [0, 1, 2]
+        assert result.nterm == 2
+
+    def test_outer_gives_union(self, m3: Model) -> None:
+        m = Model()
+        a = m.add_variables(coords=[pd.Index([0, 1, 2], name="i")], name="a")
+        b = m.add_variables(coords=[pd.Index([5, 6, 7], name="i")], name="b")
+        result = a.add(b, join="outer")
+        assert set(result.coords["i"].values) == {0, 1, 2, 5, 6, 7}
diff --git a/test/test_optimization.py b/test/test_optimization.py
index 492d703a..8f8ec6d4 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -586,13 +586,12 @@ def test_non_aligned_variables(
 
     tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL
 
+    # With inner join default, x + y covers only the intersection (dim_0=0..7).
+    # x[8] and x[9] don't appear in any constraint or objective.
     with pytest.warns(UserWarning):
         assert np.isclose(
             model_with_non_aligned_variables.solution["x"][0], 0, rtol=tol
         )
-        assert np.isclose(
-            model_with_non_aligned_variables.solution["x"][-1], 10.5, rtol=tol
-        )
         assert np.isclose(
             model_with_non_aligned_variables.solution["y"][0], 10.5, rtol=tol
         )