fixed mcwf/lindblad ordering not matching tjm (#481)

## Description

This PR fixes a problem where the MCWF/Lindblad solvers had a different
qubit ordering due to the Kronecker product. This has been fixed so that
the qubit order remains consistent for all methods across the whole
library.

## Checklist

<!---
This checklist serves as a reminder of a couple of things that ensure
your pull request will be merged swiftly.
-->

- [x] The pull request only contains commits that are focused and
relevant to this change.
- [x] I have added appropriate tests that cover the new/changed
functionality.
- [x] I have updated the documentation to reflect these changes.
- [x] I have added entries to the changelog for any noteworthy
additions, changes, fixes, or removals.
- [x] I have added migration instructions to the upgrade guide (if
needed).
- [x] The changes follow the project's style guidelines and introduce no
new warnings.
- [x] The changes are fully tested and pass the CI checks.
- [x] I have reviewed my own code changes.

**If PR contains AI-assisted content:**

- [x] I have disclosed the use of AI tools in the PR description as per
our [AI Usage
Guidelines](https://github.com/munich-quantum-toolkit/yaqs/blob/main/docs/ai_usage.md).
- [x] AI-assisted commits include an `Assisted-by: [Model Name] via
[Tool Name]` footer.
- [x] I confirm that I have personally reviewed and understood all
AI-generated content, and accept full responsibility for it.

Author: aaronleesander

CHANGELOG.md

@@ -48,6 +48,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel
 
 ### Fixed
 
+- fixed mcwf/lindblad ordering not matching tjm ([#481]) ([**@aaronleesander**])
 - added regression tests for tjm and mcwf jump probabilities ([#479]) ([**@aaronleesander**])
 - fixed NoiseModel factor order when two-site indices are normalized ([#396]) ([**@aleramos119**])
 - refactored observable handling in simulator to preserve user order ([#447]) ([**@aaronleesander**])
@@ -144,6 +145,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool
 
 <!-- PR links -->
 
+[#481]: https://github.com/munich-quantum-toolkit/yaqs/pull/481
 [#476]: https://github.com/munich-quantum-toolkit/yaqs/pull/476
 [#479]: https://github.com/munich-quantum-toolkit/yaqs/pull/479
 [#477]: https://github.com/munich-quantum-toolkit/yaqs/pull/477

UPGRADING.md

@@ -102,6 +102,37 @@ For MPS-backed analog and strong-digital runs, `result.runtime_cost`, `result.ma
 `result.total_bond` are filled automatically (aligned with `result.times` or the strong-sim layer
 grid). MCWF, Lindblad, and weak digital runs leave these as `None`.
 
+### MCWF / Lindblad operator ordering (dense backends)
+
+MCWF (`State(..., representation="vector")`) and Lindblad (`representation="density_matrix"`)
+embed jump operators and observables on the full Hilbert space using the same **site-0 LSB**
+convention as MPS `to_vec`, Qiskit little-endian circuits, and the TJM (MPO) dissipation path.
+Before this release, those dense embeddings used a different Kronecker-product order, so jump
+probabilities, observables, and cross-solver comparisons could disagree with TJM even when the
+`NoiseModel` definition looked identical.
+
+**What changed:** `_embed_operator_sparse` / `_embed_observable_sparse` (and their dense
+counterparts) now delegate to `state_utils.embed_*` helpers instead of building
+`left ⊗ op ⊗ right` with reversed tensor-leg order.
+
+**Why it matters:** MCWF, Lindblad, and TJM now agree on how a local operator on `sites=[i]` or
+adjacent `sites=[i, i+1]` is placed in the full space. Regression tests compare TJM dissipative
+norm loss to MCWF jump probabilities under lowering noise.
+
+**What you need to do:**
+
+- If you only pass standard `NoiseModel` processes (`sites`, built-in names, or matrices authored
+  for the listed site order), **no change is required**—results may shift slightly because the
+  previous ordering was incorrect.
+- If you hand-built full-space jump operators or compared MCWF/Lindblad outputs to TJM using
+  custom dense embeddings, rebuild those operators with
+  `mqt.yaqs.core.data_structures.state_utils.embed_one_site_operator`,
+  `embed_adjacent_two_site_operator`, or `embed_two_site_factors`, or pass the same local matrices
+  through `NoiseModel` and let the solvers embed them.
+- For adjacent two-site **matrix** processes, list sites in ascending order `[i, i+1]` with the
+  local matrix written for that pair order. If you pass reversed sites `[i+1, i]`, the matrix is
+  transposed automatically to match the `(i, i+1)` leg order.
+
 ### Top-level public API
 
 ```python

src/mqt/yaqs/analog/lindblad.py

@@ -112,6 +112,7 @@ def preprocess_lindblad(
         raise ValueError(msg)
 
     dim = math.prod(resolve_physical_dimensions(num_sites, physical_dimensions))
+    site_dims = resolve_physical_dimensions(num_sites, physical_dimensions)
 
     if dim > 2**10:
         msg = (
@@ -166,7 +167,7 @@ def preprocess_lindblad(
             strength = process["strength"]
             if strength <= 0:
                 continue
-            op_full = _embed_operator_sparse(process, num_sites)
+            op_full = _embed_operator_sparse(process, num_sites, physical_dimensions=site_dims)
             jump_ops.append(np.sqrt(strength) * op_full)
 
     is_unitary = len(jump_ops) == 0
@@ -184,7 +185,7 @@ def preprocess_lindblad(
         if obs.gate.name in {"entropy", "schmidt_spectrum"}:
             embedded_observables.append(None)
         else:
-            embedded_observables.append(_embed_observable_sparse(obs, num_sites))
+            embedded_observables.append(_embed_observable_sparse(obs, num_sites, physical_dimensions=site_dims))
 
     # 6. Fixed-step propagator exp(L dt) when vec(rho) fits in memory (time-independent generator).
     step_propagator: NDArray[np.complex128] | None = None

src/mqt/yaqs/analog/mcwf.py

@@ -117,6 +117,7 @@ def preprocess_mcwf(
         raise ValueError(msg)
 
     dim = math.prod(resolve_physical_dimensions(num_sites, physical_dimensions))
+    site_dims = resolve_physical_dimensions(num_sites, physical_dimensions)
 
     if dim > 2**14:
         msg = (
@@ -161,7 +162,7 @@ def preprocess_mcwf(
             strength = process["strength"]
             if strength <= 0:
                 continue
-            op_full = _embed_operator_sparse(process, num_sites)
+            op_full = _embed_operator_sparse(process, num_sites, physical_dimensions=site_dims)
             jump_ops.append(np.sqrt(strength) * op_full)
 
     is_unitary = len(jump_ops) == 0
@@ -190,7 +191,7 @@ def preprocess_mcwf(
         if obs.gate.name in {"entropy", "schmidt_spectrum"}:
             embedded_observables.append(None)
         else:
-            op = _embed_observable_sparse(obs, num_sites)
+            op = _embed_observable_sparse(obs, num_sites, physical_dimensions=site_dims)
             embedded_observables.append(op)
 
     return MCWFContext(