Add embed_lazy for lazy operator embedding

[Temitope15]

Summary

Adds a lazy operator-embedding API, embed_lazy, as requested in the unitaryHACK bounty issue qojulia/QuantumOptics.jl#523.

Where embed materializes the full embedded matrix (e.g. a 2^n × 2^n sparse matrix to embed a one-site operator into an n-site basis), embed_lazy preserves the local tensor-product structure and returns a structure-preserving lazy operator. The result can be multiplied by kets, used as a term in a LazySum, and converted with dense/sparse only when explicitly requested.

b0 = SpinBasis(1//2)
b = b0 ⊗ b0 ⊗ b0 ⊗ b0          # composite basis
H = embed_lazy(b, 1, sigmax(b0)) + embed_lazy(b, 2, sigmaz(b0))

What it does

• DataOperator / any single-subsystem AbstractOperator at a single index → LazyTensor
• LazyTensor re-embedded into the larger basis, preserving suboperators and factor (the suboperator indices are remapped through indices, so partial LazyTensors embed correctly)
• LazySum preserved as a LazySum of lazily embedded terms, factors carried over
• TimeDependentSum preserved, keeping its coefficients while lazily embedding the static terms
• single- and multi-index forms, plus the explicit embed_lazy(basis_l, basis_r, indices, op) signature
• incompatible bases raise IncompatibleBases; unsorted indices raise ArgumentError

API decisions (the points raised in the issue)

• Name: kept embed_lazy from PR Draft: Embedding with lazy semantics #88, exported from QuantumOpticsBase.
• Two-basis form: embed_lazy(basis_l, basis_r, indices, op) is supported alongside the single-basis form.
• Monolithic multi-site operators: a single dense/sparse operator acting jointly on more than one subsystem has no tensor-product structure to keep, so it cannot be embedded lazily — these raise an informative ArgumentError pointing the user to embed or to passing a LazyTensor/LazySum. (LazyProduct is likewise left unsupported for now.)

Note on the issue's example

The issue's minimal example writes local_h = sx + 0.5*sz and expects embed_lazy(b, 3, local_h) isa LazySum. In QuantumOpticsBase, + on two operators sharing a basis returns a single eager Operator, not a LazySum, so embed_lazy correctly returns a LazyTensor for that input. embed_lazy preserves a LazySum only when actually given one (e.g. LazySum(sx, 0.5*sz)); the tests and docstring use an explicit LazySum to reflect this.

Tests & docs

• New test/test_embed_lazy.jl: single-index dense/sparse, the two-basis signature, LazySum preservation + factor carry-over, multi-index LazyTensor re-embedding (including a partial LazyTensor and sparse suboperators), TimeDependentSum, incompatible-basis and unsupported-operator failures, and an allocation check confirming the lazy path does not materialize the full matrix.
• Docstring with a runnable jldoctest and an @docs entry in docs/src/api.md next to embed.
• Verified locally on Julia 1.11.5: the new test_embed_lazy test set, the Aqua checks (no new ambiguities), the doctests, and JET on the new methods all pass.

AI disclosure

I used Claude Code as an assistant to help me find my way around this large codebase and to help brainstorm and draft parts of the implementation, tests, and docstring. I drove the design decisions, reviewed and refined the code, and manually verified everything on Julia 1.11.5 — running the new test_embed_lazy set, the Aqua checks, the doctests, and JET. The clarification about + returning an eager operator (so embed_lazy returns a LazyTensor, not a LazySum, for the issue's example) and the choice to fail fast on unsorted indices both came from problems I hit and worked through while testing locally.

[Krastanov]

we do not accept LLM-generated PRs for bounties