Calling to_matrix on operator seems to be very slow

[SimonRicher]

Required prerequisites

• Consult the security policy. If reporting a security vulnerability, do not report the bug using this form. Use the process described in the policy to report the issue.
• Make sure you've read the documentation. Your issue may be addressed there.
• Search the issue tracker to verify that this hasn't already been reported. +1 or comment there if it has.
• If possible, make a PR with a failing test to give us a starting point to work on!

Describe the bug

It seems to be very slow to create the full dense matrix from a cudaQ lazy representation. It takes 13 seconds for a chain of 9 transmon qubits and 346 seconds for a chain of 10 transmon qubits on my computer.

Steps to reproduce the bug

import cudaq
from cudaq import operators
import numpy as np
import timeit 
import scqubits as scq

cudaq.set_target('dynamics')
GHz = 2 * np.pi
MHz = 2 * np.pi * 1e-3


def create_transmon(Ej: float, Ec: float,dimension:int)->tuple[np.array,np.array]:

    ### Create transmon in fock basis
    transmon = scq.Transmon(
     # Number of energy levels to include
    EC = Ec * GHz,        # Charging energy (Hz)
    EJ = Ej * GHz,        # Josephson energy of junction (Hz)
    ng = 0,
    ncut = 1000,
    truncated_dim=dimension,
    )
    return transmon.hamiltonian(energy_esys=True).astype(np.complex128),transmon.n_operator(energy_esys=True).astype(np.complex128)

chain_len, qubit_dim = 10,2

Ejs = 20
Ecs = 0.1
dimensions = dict(enumerate([qubit_dim]*chain_len))

transmon_op = [create_transmon(Ejs,Ecs,qubit_dim) for _ in range(chain_len)]
n = []
H = []

for i in range(chain_len):
    H_i, n_i = transmon_op[i]
    operators.define(f"H_{i}", [qubit_dim], lambda H=H_i:H)
    H.append(operators.instantiate(f"H_{i}",[i]))
    operators.define(f"n_{i}", [qubit_dim], lambda n=n_i : n)
    n.append(operators.instantiate(f"n_{i}",[i]))
H_chain = H[0]
for i in range(chain_len):
    H_chain += H[i]
    if i < chain_len-1:
        H_chain +=  n[i] * n[i+1]

### Benchmark the matrix materialization
number = 1
code_to_benchmark = """H_chain.to_matrix(dimensions)
"""
runtime = timeit.timeit(code_to_benchmark,globals = globals(), number = number)
print(f"Average time cudaq : {runtime/number}",flush = True)

Expected behavior

I don't think it should be that long to go from a lazy representation to a dense representation for this size of matrix (1024X1024).

Is this a regression? If it is, put the last known working version (or commit) here.

Not a regression

Environment

• CUDA-Q version: 0.13.0
• Python version: 3.12.9
• C++ compiler:
• Operating system: Ubuntu 22.04.5 LTS

Suggestions

No response