Base Gaussian State dm doesn't preserve locality dimensions

Before posting a bug report

[X] I have searched exisisting GitHub issues to make sure the issue does not already exist.

Expected behavior

BaseGaussianState.dm() doesn't preserve locality dimensions. The Fock basis backend and Gaussian backend should give the same results.

Actual behavior

state.dm() gives gives only 2 dimensional matricies. The Fock basis state gives a multi-index object depending on the number of modes.

Reproduces how often

Always

System information

>>> sf.about()

Strawberry Fields: a Python library for continuous-variable quantum circuits.
Copyright 2018-2020 Xanadu Quantum Technologies Inc.

Python version:            3.10.5
Platform info:             Linux-5.18.12-arch1-1-x86_64-with-glibc2.35
Installation path:         /home/lund/.local/lib/python3.10/site-packages/strawberryfields
Strawberry Fields version: 0.23.0
Numpy version:             1.22.4
Scipy version:             1.8.1
SymPy version:             1.10.1
NetworkX version:          2.8.5
The Walrus version:        0.19.0
Blackbird version:         0.4.0
XCC version:               0.2.1
TensorFlow version:        None

Source code

import strawberryfields as sf
import strawberryfields.ops as op

p = sf.Program(2)
with p.context as q:
    op.Squeezed(1) | q[0]

e = sf.Engine('fock', backend_options={'cutoff_dim': 4})
r = e.run(p)
print(f"Fock engine density matrix shape {r.state.dm().shape}")

e = sf.Engine('gaussian')
r = e.run(p)
print(f"Gaussian engine density matrix shape {r.state.dm(cutoff=4).shape}")

Tracebacks

No response

Additional information