Repository: PennyLaneAI/pennylane (python)
Issue: QNodes with `default.qubit` and `default.qubit.legacy` give different results 
### Expected behavior

They should give the same results.

### Actual behavior

They give different results.

### Additional information

Seems to be that BasisState followed by some (? — RandomLayers and some instances of QubitUnitary ) operation give different results.

### Source code

```shell
import pennylane as qml
import pennylane.numpy as np

new_device = qml.device("default.qubit", wires=2)
legacy_device = qml.device("default.qubit.legacy", wires=2)

def random_unitary():
    weights = np.array([[0.1, -2.1, 1.4]])
    qml.RandomLayers(weights=weights, wires=0)

def circuit(state):
    qml.BasisState(state, wires=1) 
    random_unitary()
    # getting rid of BasisState and/or random_unitary produces no discrepancy

@qml.qnode(new_device)
def circuit_new_device(state):
    circuit(state)
    return qml.state()

@qml.qnode(legacy_device)
def circuit_legacy_device(state):
    circuit(state)
    return qml.state()

res_new_device = circuit_new_device(np.array([1]))
res_legacy_device = circuit_legacy_device(np.array([1]))

print(res_new_device == res_legacy_device) # [False False False False]
```


### Tracebacks

```shell
N/A
```


### System information

```shell
Name: PennyLane
Version: 0.33.0
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/PennyLaneAI/pennylane
Author:
Author-email:
License: Apache License 2.0
Location: /Users/isaac/.virtualenvs/pennylane-stable/lib/python3.9/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane-Lightning

Platform info:           macOS-14.1-x86_64-i386-64bit
Python version:          3.9.14
Numpy version:           1.23.5
Scipy version:           1.10.0
Installed devices:
- default.gaussian (PennyLane-0.34.0.dev0)
- default.mixed (PennyLane-0.34.0.dev0)
- default.qubit (PennyLane-0.34.0.dev0)
- default.qubit.autograd (PennyLane-0.34.0.dev0)
- default.qubit.jax (PennyLane-0.34.0.dev0)
- default.qubit.legacy (PennyLane-0.34.0.dev0)
- default.qubit.tf (PennyLane-0.34.0.dev0)
- default.qubit.torch (PennyLane-0.34.0.dev0)
- default.qutrit (PennyLane-0.34.0.dev0)
- null.qubit (PennyLane-0.34.0.dev0)
- lightning.qubit (PennyLane-Lightning-0.33.1)
```


### Existing GitHub issues

- [X] I have searched existing GitHub issues to make sure the issue does not already exist.
Pull request description: Fix the wire order of state prep operations in `default.qubit`
Fixes #4779 . [sc-49381]

When creating an initial state, `get_final_state` was using `tape.op_wires` for the wire order, even though the simulation always occurred using the wire order `(0,1,2,3,...)`.  Sometimes `tape.op_wires` would have a different order, leading to a mismatch.

Now we use `sorted(circuit.op_wires)`, so they will always be in the monotonically increasing integers used by the simulation.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 889b723d4..e9a055ae3 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -51,6 +51,10 @@
   is a `Hamiltonian`.
   [(#4768)](https://github.com/PennyLaneAI/pennylane/pull/4768)
 
+* In `default.qubit`, initial states are now initialized with the simulator's wire order, not the circuit's
+  wire order.
+  [(#4781)](https://github.com/PennyLaneAI/pennylane/pull/4781)
+
 <h3>Contributors ✍️</h3>
 
 This release contains contributions from (in alphabetical order):
diff --git a/pennylane/devices/qubit/simulate.py b/pennylane/devices/qubit/simulate.py
index 067a8fe22..22b67e474 100644
--- a/pennylane/devices/qubit/simulate.py
+++ b/pennylane/devices/qubit/simulate.py
@@ -153,7 +153,7 @@ def get_final_state(circuit, debugger=None, interface=None):
     if len(circuit) > 0 and isinstance(circuit[0], qml.operation.StatePrepBase):
         prep = circuit[0]
 
-    state = create_initial_state(circuit.op_wires, prep, like=INTERFACE_TO_LIKE[interface])
+    state = create_initial_state(sorted(circuit.op_wires), prep, like=INTERFACE_TO_LIKE[interface])
 
     # initial state is batched only if the state preparation (if it exists) is batched
     is_state_batched = bool(prep and prep.batch_size is not None)


Test patch:
diff --git a/tests/devices/default_qubit/test_default_qubit.py b/tests/devices/default_qubit/test_default_qubit.py
index 257d57089..a7178c77a 100644
--- a/tests/devices/default_qubit/test_default_qubit.py
+++ b/tests/devices/default_qubit/test_default_qubit.py
@@ -317,6 +317,17 @@ class TestBasicCircuit:
         assert qml.math.get_interface(res) == "tensorflow"
         assert qml.math.allclose(res, -1)
 
+    def test_basis_state_wire_order(self):
+        """Test that the wire order is correct with a basis state if the tape wires have a non standard order."""
+
+        dev = DefaultQubit()
+
+        tape = qml.tape.QuantumScript([qml.BasisState([1], wires=1), qml.PauliZ(0)], [qml.state()])
+
+        expected = np.array([0, 1, 0, 0], dtype=np.complex128)
+        res = dev.execute(tape)
+        assert qml.math.allclose(res, expected)
+
 
 class TestSampleMeasurements:
     """A copy of the `qubit.simulate` tests, but using the device"""