diff --git a/qiskit/quantum_info/operators/symplectic/pauli_list.py b/qiskit/quantum_info/operators/symplectic/pauli_list.py
index 02f9a78052b7..7210a44a7bbb 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli_list.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli_list.py
@@ -445,16 +445,14 @@ def insert(self, ind: int, value: PauliList, qubit: bool = False) -> PauliList:
                 f"Index {ind} is greater than number of qubits"
                 f" in the PauliList ({self.num_qubits})"
             )
-        if len(value) == 1:
-            # Pad blocks to correct size
-            value_x = np.vstack(size * [value.x])
-            value_z = np.vstack(size * [value.z])
-            value_phase = np.vstack(size * [value.phase])
-        elif len(value) == size:
+        if len(value) == size:
             #  Blocks are already correct size
             value_x = value.x
             value_z = value.z
-            value_phase = value.phase
+        elif len(value) == 1:
+            # Pad blocks to correct size
+            value_x = np.vstack(size * [value.x])
+            value_z = np.vstack(size * [value.z])
         else:
             # Blocks are incorrect size
             raise QiskitError(
@@ -465,7 +463,7 @@ def insert(self, ind: int, value: PauliList, qubit: bool = False) -> PauliList:
         # Build new array by blocks
         z = np.hstack([self.z[:, :ind], value_z, self.z[:, ind:]])
         x = np.hstack([self.x[:, :ind], value_x, self.x[:, ind:]])
-        phase = self.phase + value_phase
+        phase = self.phase + value.phase
 
         return PauliList.from_symplectic(z, x, phase)
 
@@ -1121,6 +1119,8 @@ def from_symplectic(
         Returns:
             PauliList: the constructed PauliList.
         """
+        if isinstance(phase, np.ndarray) and np.ndim(phase) > 1:
+            raise ValueError(f"phase should be at most 1D but has {np.ndim(phase)} dimensions.")
         base_z, base_x, base_phase = cls._from_array(z, x, phase)
         return cls(BasePauli(base_z, base_x, base_phase))
 
diff --git a/releasenotes/notes/fix-paulilist-length1-phase-688d0e3a64ec9a9f.yaml b/releasenotes/notes/fix-paulilist-length1-phase-688d0e3a64ec9a9f.yaml
new file mode 100644
index 000000000000..a1a59708003c
--- /dev/null
+++ b/releasenotes/notes/fix-paulilist-length1-phase-688d0e3a64ec9a9f.yaml
@@ -0,0 +1,6 @@
+---
+fixes:
+  - |
+    Fixed a bug that caused :meth:`.PauliList.insert` with ``qubit=True`` to produce a `phase`
+    attribute with the wrong shape when the original object was length 1.
+    Fixed `#13623 <https://github.com/Qiskit/qiskit/issues/13623>`__.
diff --git a/test/python/quantum_info/operators/symplectic/test_pauli_list.py b/test/python/quantum_info/operators/symplectic/test_pauli_list.py
index 8c96f63c4ddf..9abc473dc333 100644
--- a/test/python/quantum_info/operators/symplectic/test_pauli_list.py
+++ b/test/python/quantum_info/operators/symplectic/test_pauli_list.py
@@ -1560,6 +1560,15 @@ def test_insert(self):
                 value1 = pauli.insert(1, insert)
                 self.assertEqual(value1, target1)
 
+        # Insert single column to length-1 PauliList:
+        with self.subTest(msg="length-1, single-column, single-val"):
+            pauli = PauliList(["X"])
+            insert = PauliList(["Y"])
+            target0 = PauliList(["YX"])
+            value0 = pauli.insert(1, insert, qubit=True)
+            self.assertEqual(value0, target0)
+            self.assertEqual(value0.phase.shape, (1,))
+
         # Insert single column
         pauli = PauliList(["X", "Y", "Z", "-iI"])
         for i in ["I", "X", "Y", "Z", "iY"]: