diff --git a/compiler/qsc/src/bin/qsi.rs b/compiler/qsc/src/bin/qsi.rs
index b316746364..d2f7a455e6 100644
--- a/compiler/qsc/src/bin/qsi.rs
+++ b/compiler/qsc/src/bin/qsi.rs
@@ -81,6 +81,16 @@ impl Receiver for TerminalReceiver {
         Ok(())
     }
 
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> std::result::Result<(), output::Error> {
+        println!("Matrix:");
+        for row in matrix {
+            let row = row.iter().map(|elem| format!("[{}, {}]", elem.re, elem.im));
+            println!("{}", row.collect::<Vec<_>>().join(", "));
+        }
+
+        Ok(())
+    }
+
     fn message(&mut self, msg: &str) -> Result<(), output::Error> {
         println!("{msg}");
         Ok(())
diff --git a/compiler/qsc/src/lib.rs b/compiler/qsc/src/lib.rs
index 65029ed79a..db9ee9ff4c 100644
--- a/compiler/qsc/src/lib.rs
+++ b/compiler/qsc/src/lib.rs
@@ -49,7 +49,10 @@ pub mod line_column {
 
 pub use qsc_eval::{
     backend::{Backend, SparseSim},
-    state::{fmt_basis_state_label, fmt_complex, format_state_id, get_latex, get_phase},
+    state::{
+        fmt_basis_state_label, fmt_complex, format_state_id, get_matrix_latex, get_phase,
+        get_state_latex,
+    },
 };
 
 pub mod linter {
diff --git a/compiler/qsc_eval/src/intrinsic.rs b/compiler/qsc_eval/src/intrinsic.rs
index 24498d970f..274df93578 100644
--- a/compiler/qsc_eval/src/intrinsic.rs
+++ b/compiler/qsc_eval/src/intrinsic.rs
@@ -77,6 +77,24 @@ pub(crate) fn call(
                 Err(_) => Err(Error::OutputFail(name_span)),
             }
         }
+        "DumpMatrix" => {
+            let qubits = arg.unwrap_array();
+            let qubits = qubits
+                .iter()
+                .map(|q| q.clone().unwrap_qubit().0)
+                .collect::<Vec<_>>();
+            if qubits.len() != qubits.iter().collect::<FxHashSet<_>>().len() {
+                return Err(Error::QubitUniqueness(arg_span));
+            }
+            let (state, qubit_count) = sim.capture_quantum_state();
+            let state = utils::split_state(&qubits, &state, qubit_count)
+                .map_err(|()| Error::QubitsNotSeparable(arg_span))?;
+            let matrix = utils::state_to_matrix(state, qubits.len() / 2);
+            match out.matrix(matrix) {
+                Ok(()) => Ok(Value::unit()),
+                Err(_) => Err(Error::OutputFail(name_span)),
+            }
+        }
         "PermuteLabels" => qubit_relabel(arg, arg_span, |q0, q1| sim.qubit_swap_id(q0, q1)),
         "Message" => match out.message(&arg.unwrap_string()) {
             Ok(()) => Ok(Value::unit()),
diff --git a/compiler/qsc_eval/src/intrinsic/tests.rs b/compiler/qsc_eval/src/intrinsic/tests.rs
index d06d35e8b4..34b4291480 100644
--- a/compiler/qsc_eval/src/intrinsic/tests.rs
+++ b/compiler/qsc_eval/src/intrinsic/tests.rs
@@ -440,6 +440,26 @@ fn dump_register_other_qubits_one_state_is_separable() {
     );
 }
 
+#[test]
+fn dump_register_other_qubits_phase_reflected_in_subset() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[3];
+            H(qs[0]);
+            X(qs[2]);
+            Z(qs[2]);
+            Microsoft.Quantum.Diagnostics.DumpRegister(qs[...1]);
+            ResetAll(qs);
+        }"},
+        &expect![[r#"
+            STATE:
+            |00⟩: −0.7071+0.0000𝑖
+            |10⟩: −0.7071+0.0000𝑖
+        "#]],
+    );
+}
+
 #[test]
 fn dump_register_qubits_reorder_output() {
     check_intrinsic_output(
diff --git a/compiler/qsc_eval/src/intrinsic/utils.rs b/compiler/qsc_eval/src/intrinsic/utils.rs
index 79020689da..ed1f38c4d9 100644
--- a/compiler/qsc_eval/src/intrinsic/utils.rs
+++ b/compiler/qsc_eval/src/intrinsic/utils.rs
@@ -5,8 +5,8 @@ use std::collections::hash_map::Entry;
 
 use num_bigint::BigUint;
 use num_complex::{Complex, Complex64};
-use num_traits::{One, Zero};
-use rustc_hash::FxHashMap;
+use num_traits::Zero;
+use rustc_hash::{FxHashMap, FxHashSet};
 
 /// Given a state and a set of qubits, split the state into two parts: the qubits to dump and the remaining qubits.
 /// This function will return an error if the state is not separable using the provided qubit identifiers.
@@ -22,26 +22,13 @@ pub fn split_state(
     }
 
     let mut dump_state = FxHashMap::default();
-    let mut other_state = FxHashMap::default();
 
     // Compute the mask for the qubits to dump and the mask for the other qubits.
     let (dump_mask, other_mask) = compute_mask(qubit_count, qubits);
 
     // Try to split out the state for the given qubits from the whole state, detecting any entanglement
     // and returning an error if the qubits are not separable.
-    let dump_norm = collect_split_state(
-        state,
-        &dump_mask,
-        &other_mask,
-        &mut dump_state,
-        &mut other_state,
-    )?;
-
-    // If the product of the collected states is not equal to the total number of input states, then that
-    // implies some states are zero amplitude that would have to be non-zero for the state to be separable.
-    if state.len() != dump_state.len() * other_state.len() {
-        return Err(());
-    }
+    let dump_norm = collect_split_state(state, &dump_mask, &other_mask, &mut dump_state)?;
 
     let dump_norm = 1.0 / dump_norm.sqrt();
     let mut dump_state = dump_state
@@ -79,7 +66,6 @@ fn collect_split_state(
     dump_mask: &BigUint,
     other_mask: &BigUint,
     dump_state: &mut FxHashMap<BigUint, Complex64>,
-    other_state: &mut FxHashMap<BigUint, Complex64>,
 ) -> Result<f64, ()> {
     // To ensure consistent ordering, we iterate over the vector directly (returned from the simulator in a deterministic order),
     // and not the map used for arbitrary lookup.
@@ -88,12 +74,11 @@ fn collect_split_state(
     let (base_label, base_val) = state_iter.next().expect("state should never be empty");
     let dump_base_label = base_label & dump_mask;
     let other_base_label = base_label & other_mask;
-    let mut dump_norm = 1.0_f64;
+    let mut dump_norm = base_val.norm().powi(2);
+    let mut other_state = FxHashSet::default();
 
-    // Start with an amplitude of 1 in the first expected split states. This becomes the basis
-    // of the split later, but will get normalized as part of collecting the remaining states.
-    dump_state.insert(dump_base_label.clone(), Complex64::one());
-    other_state.insert(other_base_label.clone(), Complex64::one());
+    dump_state.insert(dump_base_label.clone(), *base_val);
+    other_state.insert(other_base_label.clone());
 
     for (curr_label, curr_val) in state_iter {
         let dump_label = curr_label & dump_mask;
@@ -117,23 +102,23 @@ fn collect_split_state(
         }
 
         if let Entry::Vacant(entry) = dump_state.entry(dump_label) {
-            // When capturing the amplitude for the dump state, we must divide out the amplitude for the other
-            // state, and vice-versa below.
-            let amplitude = curr_val / other_val;
+            let amplitude = *curr_val;
             let norm = amplitude.norm().powi(2);
             if !norm.is_nearly_zero() {
                 entry.insert(amplitude);
                 dump_norm += norm;
             }
         }
-        if let Entry::Vacant(entry) = other_state.entry(other_label) {
-            let amplitude = curr_val / dump_val;
-            let norm = amplitude.norm().powi(2);
-            if !norm.is_nearly_zero() {
-                entry.insert(amplitude);
-            }
+        if !(curr_val / dump_val).norm().powi(2).is_nearly_zero() {
+            other_state.insert(other_label);
         }
     }
+
+    // If the product of the collected states is not equal to the total number of input states, then that
+    // implies some states are zero amplitude that would have to be non-zero for the state to be separable.
+    if state.len() != dump_state.len() * other_state.len() {
+        return Err(());
+    }
     Ok(dump_norm)
 }
 
@@ -185,3 +170,28 @@ where
         self.re.is_nearly_zero() && self.im.is_nearly_zero()
     }
 }
+
+pub(crate) fn state_to_matrix(
+    state: Vec<(BigUint, Complex64)>,
+    qubit_count: usize,
+) -> Vec<Vec<Complex64>> {
+    let state: FxHashMap<BigUint, Complex<f64>> = state.into_iter().collect();
+    let mut matrix = Vec::new();
+    let num_entries: usize = 1 << qubit_count;
+    #[allow(clippy::cast_precision_loss)]
+    let factor = (num_entries as f64).sqrt();
+    for i in 0..num_entries {
+        let mut row = Vec::new();
+        for j in 0..num_entries {
+            let key = BigUint::from(i * num_entries + j);
+            let val = match state.get(&key) {
+                Some(val) => val * factor,
+                None => Complex::zero(),
+            };
+            row.push(val);
+        }
+        matrix.push(row);
+    }
+
+    matrix
+}
diff --git a/compiler/qsc_eval/src/output.rs b/compiler/qsc_eval/src/output.rs
index 4c4f569f34..cb4a9cced7 100644
--- a/compiler/qsc_eval/src/output.rs
+++ b/compiler/qsc_eval/src/output.rs
@@ -16,6 +16,11 @@ pub trait Receiver {
     /// This will return an error if handling the output fails.
     fn state(&mut self, state: Vec<(BigUint, Complex64)>, qubit_count: usize) -> Result<(), Error>;
 
+    /// Receive matrix output
+    /// # Errors
+    /// This will return an error if handling the output fails.
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> Result<(), Error>;
+
     /// Receive generic message output
     /// # Errors
     /// This will return an error if handling the output fails.
@@ -47,6 +52,15 @@ impl<'a> Receiver for GenericReceiver<'a> {
         Ok(())
     }
 
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> Result<(), Error> {
+        writeln!(self.writer, "MATRIX:").map_err(|_| Error)?;
+        for row in matrix {
+            let row_str = row.iter().map(fmt_complex).collect::<Vec<_>>().join(" ");
+            writeln!(self.writer, "{row_str}").map_err(|_| Error)?;
+        }
+        Ok(())
+    }
+
     fn message(&mut self, msg: &str) -> Result<(), Error> {
         writeln!(self.writer, "{msg}").map_err(|_| Error)
     }
@@ -86,6 +100,15 @@ impl<'a> Receiver for CursorReceiver<'a> {
         Ok(())
     }
 
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> Result<(), Error> {
+        writeln!(self.cursor, "MATRIX:").map_err(|_| Error)?;
+        for row in matrix {
+            let row_str = row.iter().map(fmt_complex).collect::<Vec<_>>().join(" ");
+            writeln!(self.cursor, "{row_str}").map_err(|_| Error)?;
+        }
+        Ok(())
+    }
+
     fn message(&mut self, msg: &str) -> Result<(), Error> {
         writeln!(self.cursor, "{msg}").map_err(|_| Error)
     }
diff --git a/compiler/qsc_eval/src/state.rs b/compiler/qsc_eval/src/state.rs
index 6e69c071b3..85dc2f524c 100644
--- a/compiler/qsc_eval/src/state.rs
+++ b/compiler/qsc_eval/src/state.rs
@@ -6,7 +6,7 @@ mod tests;
 
 use num_bigint::BigUint;
 use num_complex::{Complex, Complex64};
-use std::fmt::Write;
+use std::{f64::consts::FRAC_1_SQRT_2, fmt::Write};
 
 #[must_use]
 pub fn format_state_id(id: &BigUint, qubit_count: usize) -> String {
@@ -209,7 +209,7 @@ impl RealNumber {
     }
 }
 
-/// Represents a non-zero complex numbers in the polar form: ``coefficient·𝒆^(𝝅·𝒊·phase_multiplier)``
+/// Represents a non-zero complex numbers in the polar form: ``magnitude·𝒆^(𝝅·𝒊·phase_multiplier)``
 /// Sign of the number is separated for easier composition and rendering
 #[derive(Debug)]
 struct PolarForm {
@@ -248,6 +248,10 @@ impl PolarForm {
 
     /// Try to recognize a complex number and represent it in the polar form.
     fn recognize(re: f64, im: f64) -> Option<Self> {
+        if !is_significant(re.abs()) && !is_significant(im.abs()) {
+            // 0 is better represented in Cartesian form, not polar.
+            return None;
+        }
         for (pi_num, pi_den) in Self::PI_FRACTIONS {
             #[allow(clippy::cast_precision_loss)] // We only use fixes set of fractions
             let angle: f64 = std::f64::consts::PI * (pi_num as f64) / (pi_den as f64);
@@ -345,7 +349,7 @@ fn get_terms_for_state(state: &Vec<(BigUint, Complex64)>) -> Vec<Term> {
 /// `None` is returned if the resulting formula is not nice, i.e.
 /// if the formula consists of more than 16 terms or if more than two coefficients are not recognized.
 #[must_use]
-pub fn get_latex(state: &Vec<(BigUint, Complex64)>, qubit_count: usize) -> Option<String> {
+pub fn get_state_latex(state: &Vec<(BigUint, Complex64)>, qubit_count: usize) -> Option<String> {
     if state.len() > 16 {
         return None;
     }
@@ -381,13 +385,94 @@ pub fn get_latex(state: &Vec<(BigUint, Complex64)>, qubit_count: usize) -> Optio
     Some(latex)
 }
 
+fn is_close_enough(val: &Complex64, target: &Complex64) -> bool {
+    (val.re - target.re).abs() < 1e-9 && (val.im - target.im).abs() < 1e-9
+}
+
+// Quick and dirty matching for the most common matrix elements we care about rendering
+// LaTeX for, e.g.  1/sqrt(2), -i/sqrt(2), etc.
+// Anything not in this list gets a standard rendering.
+fn get_latex_for_simple_term(val: &Complex64) -> Option<String> {
+    if is_close_enough(val, &Complex64::new(FRAC_1_SQRT_2, 0.0)) {
+        return Some("\\frac{1}{\\sqrt{2}}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.0, FRAC_1_SQRT_2)) {
+        return Some("\\frac{i}{\\sqrt{2}}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(-FRAC_1_SQRT_2, 0.0)) {
+        return Some("-\\frac{1}{\\sqrt{2}}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.0, -FRAC_1_SQRT_2)) {
+        return Some("-\\frac{i}{\\sqrt{2}}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.0, 0.5)) {
+        return Some("\\frac{i}{2}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.0, -0.5)) {
+        return Some("-\\frac{i}{2}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.5, 0.5)) {
+        return Some("\\frac{1}{2} + \\frac{i}{2}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(-0.5, -0.5)) {
+        return Some("-\\frac{1}{2} - \\frac{i}{2}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(-0.5, 0.5)) {
+        return Some("-\\frac{1}{2} + \\frac{i}{2}".to_string());
+    }
+    if is_close_enough(val, &Complex64::new(0.5, -0.5)) {
+        return Some("\\frac{1}{2} - \\frac{i}{2}".to_string());
+    }
+    None
+}
+
+#[must_use]
+pub fn get_matrix_latex(matrix: &Vec<Vec<Complex64>>) -> String {
+    let mut latex: String = String::with_capacity(500);
+    latex.push_str("$ \\begin{bmatrix} ");
+    for row in matrix {
+        let mut is_first: bool = true;
+        for element in row {
+            if !is_first {
+                latex.push_str(" & ");
+            }
+            is_first = false;
+
+            if let Some(simple_latex) = get_latex_for_simple_term(element) {
+                latex.push_str(&simple_latex);
+                continue;
+            }
+
+            let cpl = ComplexNumber::recognize(element.re, element.im);
+            write_latex_for_complex_number(&mut latex, &cpl);
+        }
+        latex.push_str(" \\\\ ");
+    }
+    latex.push_str("\\end{bmatrix} $");
+    latex.shrink_to_fit();
+    latex
+}
+
+/// Write latex for a standalone complex number
+/// '-', 0 and 1 are always rendered, but '+' is not.
+fn write_latex_for_complex_number(latex: &mut String, number: &ComplexNumber) {
+    match number {
+        ComplexNumber::Cartesian(cartesian_form) => {
+            write_latex_for_cartesian_form(latex, cartesian_form, false, true);
+        }
+        ComplexNumber::Polar(polar_form) => {
+            write_latex_for_polar_form(latex, polar_form, false);
+        }
+    }
+}
+
 /// Write latex for one term of quantum state.
 /// Latex is rendered for coefficient only (not for basis vector).
 /// + is rendered only if ``render_plus`` is true.
 fn write_latex_for_term(latex: &mut String, term: &Term, render_plus: bool) {
     match &term.coordinate {
         ComplexNumber::Cartesian(cartesian_form) => {
-            write_latex_for_cartesian_form(latex, cartesian_form, render_plus);
+            write_latex_for_cartesian_form(latex, cartesian_form, render_plus, false);
         }
         ComplexNumber::Polar(polar_form) => {
             write_latex_for_polar_form(latex, polar_form, render_plus);
@@ -425,11 +510,13 @@ fn write_latex_for_polar_form(latex: &mut String, polar_form: &PolarForm, render
 /// Brackets are used if both real and imaginary parts are present.
 /// If only one part is present, its sign is used as common.
 /// If both components are present, real part sign is used as common.
-/// 1 is not rendered, but + is rendered if ``render_plus`` is true.
+/// 1 is rendered if ``render_one`` is true
+/// + is rendered if ``render_plus`` is true.
 fn write_latex_for_cartesian_form(
     latex: &mut String,
     cartesian_form: &CartesianForm,
     render_plus: bool,
+    render_one: bool,
 ) {
     if cartesian_form.sign < 0 {
         latex.push('-');
@@ -438,7 +525,6 @@ fn write_latex_for_cartesian_form(
     }
     if let RealNumber::Zero = cartesian_form.real_part {
         if let RealNumber::Zero = cartesian_form.imaginary_part {
-            // NOTE: This branch is never used.
             latex.push('0');
         } else {
             // Only imaginary part present
@@ -447,7 +533,7 @@ fn write_latex_for_cartesian_form(
         }
     } else if let RealNumber::Zero = cartesian_form.imaginary_part {
         // Only real part present
-        write_latex_for_real_number(latex, &cartesian_form.real_part, false);
+        write_latex_for_real_number(latex, &cartesian_form.real_part, render_one);
     } else {
         // Both real and imaginary parts present
         latex.push_str("\\left( ");
@@ -463,7 +549,7 @@ fn write_latex_for_cartesian_form(
 }
 
 /// Write latex for real number. Note that the sign is not rendered.
-/// 1 is only rendered if ``render_one`` is true. 0 is rendered, but not used in current code.
+/// 1 is only rendered if ``render_one`` is true.
 fn write_latex_for_real_number(latex: &mut String, number: &RealNumber, render_one: bool) {
     match number {
         RealNumber::Algebraic(algebraic_number) => {
@@ -473,7 +559,6 @@ fn write_latex_for_real_number(latex: &mut String, number: &RealNumber, render_o
             write_latex_for_decimal_number(latex, decimal_number, render_one);
         }
         RealNumber::Zero => {
-            // Note: this arm is not used.
             latex.push('0');
         }
     }
diff --git a/compiler/qsc_eval/src/state/tests.rs b/compiler/qsc_eval/src/state/tests.rs
index 1f742f74fe..bcc5005d85 100644
--- a/compiler/qsc_eval/src/state/tests.rs
+++ b/compiler/qsc_eval/src/state/tests.rs
@@ -4,15 +4,18 @@
 #![allow(clippy::needless_raw_string_hashes)]
 
 use super::{
-    get_latex, write_latex_for_algebraic_number, write_latex_for_cartesian_form,
-    write_latex_for_decimal_number, write_latex_for_polar_form, write_latex_for_real_number,
-    write_latex_for_term, AlgebraicNumber, CartesianForm, ComplexNumber, DecimalNumber, PolarForm,
-    RationalNumber, RealNumber, Term,
+    get_matrix_latex, get_state_latex, write_latex_for_algebraic_number,
+    write_latex_for_cartesian_form, write_latex_for_complex_number, write_latex_for_decimal_number,
+    write_latex_for_polar_form, write_latex_for_real_number, write_latex_for_term, AlgebraicNumber,
+    CartesianForm, ComplexNumber, DecimalNumber, PolarForm, RationalNumber, RealNumber, Term,
 };
 use crate::state::{is_fractional_part_significant, is_significant};
 use expect_test::{expect, Expect};
 use num_complex::Complex64;
-use std::{f64::consts::PI, time::Instant};
+use std::{
+    f64::consts::{FRAC_1_SQRT_2, PI},
+    time::Instant,
+};
 
 #[test]
 fn check_is_significant() {
@@ -511,7 +514,7 @@ fn check_get_latex_for_real() {
 fn assert_latex_for_cartesian(expected: &Expect, re: f64, im: f64, render_plus: bool) {
     let number = CartesianForm::recognize(re, im);
     let mut latex = String::with_capacity(50);
-    write_latex_for_cartesian_form(&mut latex, &number, render_plus);
+    write_latex_for_cartesian_form(&mut latex, &number, render_plus, false);
     expected.assert_eq(&latex);
 }
 
@@ -638,22 +641,86 @@ fn check_get_latex_for_term() {
     );
 }
 
+fn assert_latex_for_complex_number(expected: &Expect, re: f64, im: f64) {
+    let n: ComplexNumber = ComplexNumber::recognize(re, im);
+    let mut latex = String::with_capacity(50);
+    write_latex_for_complex_number(&mut latex, &n);
+    expected.assert_eq(&latex);
+}
+
+#[test]
+fn check_get_latex_for_complex_number() {
+    // Future work:
+    // While rendering is correct, a better way may be the following:
+    // -(1-i) -> -1+i remove brackets for standalone number
+    // 1/2 i -> i/2
+    // √2/2 -> 1/√2
+    assert_latex_for_complex_number(&expect!([r"0"]), 0.0, 0.0);
+
+    assert_latex_for_complex_number(&expect!([r"1"]), 1.0, 0.0);
+    assert_latex_for_complex_number(&expect!([r"-1"]), -1.0, 0.0);
+    assert_latex_for_complex_number(&expect!([r"i"]), 0.0, 1.0);
+    assert_latex_for_complex_number(&expect!([r"-i"]), 0.0, -1.0);
+
+    assert_latex_for_complex_number(&expect!([r"\frac{1}{2}"]), 0.5, 0.0);
+    assert_latex_for_complex_number(&expect!([r"-\frac{1}{2}"]), -0.5, 0.0);
+    assert_latex_for_complex_number(&expect!([r"\frac{1}{2}i"]), 0.0, 0.5);
+    assert_latex_for_complex_number(&expect!([r"-\frac{1}{2}i"]), 0.0, -0.5);
+
+    assert_latex_for_complex_number(
+        &expect!([r#"\left( \frac{1}{2}+\frac{1}{2}i \right)"#]),
+        0.5,
+        0.5,
+    );
+    assert_latex_for_complex_number(
+        &expect!([r#"-\left( \frac{1}{2}-\frac{1}{2}i \right)"#]),
+        -0.5,
+        0.5,
+    );
+    assert_latex_for_complex_number(
+        &expect!([r#"\left( \frac{1}{2}-\frac{1}{2}i \right)"#]),
+        0.5,
+        -0.5,
+    );
+    assert_latex_for_complex_number(
+        &expect!([r#"-\left( \frac{1}{2}+\frac{1}{2}i \right)"#]),
+        -0.5,
+        -0.5,
+    );
+
+    assert_latex_for_complex_number(&expect!([r#"\frac{\sqrt{2}}{2}"#]), FRAC_1_SQRT_2, 0.0);
+    assert_latex_for_complex_number(&expect!([r#"-\frac{\sqrt{2}}{2}"#]), -FRAC_1_SQRT_2, 0.0);
+    assert_latex_for_complex_number(&expect!([r#"\frac{\sqrt{2}}{2}i"#]), 0.0, FRAC_1_SQRT_2);
+    assert_latex_for_complex_number(&expect!([r#"-\frac{\sqrt{2}}{2}i"#]), 0.0, -FRAC_1_SQRT_2);
+
+    assert_latex_for_complex_number(
+        &expect!([r"\frac{1}{2} e^{ i \pi / 3}"]),
+        1.0 / 2.0 * (PI / 3.0).cos(),
+        1.0 / 2.0 * (PI / 3.0).sin(),
+    );
+    assert_latex_for_complex_number(
+        &expect!([r#"\left( \frac{1}{2}+\frac{1}{2}i \right)"#]),
+        1.0 / 2.0,
+        1.0 / 2.0,
+    );
+}
+
 #[test]
 fn check_get_latex() {
     expect!([r"$|\psi\rangle = \left( \frac{1}{2}+\frac{1}{2}i \right)|00\rangle$"]).assert_eq(
-        &get_latex(&vec![(0_u8.into(), Complex64::new(0.5, 0.5))], 2)
+        &get_state_latex(&vec![(0_u8.into(), Complex64::new(0.5, 0.5))], 2)
             .expect("expected valid latex"),
     );
     expect!([r"$|\psi\rangle = -|00\rangle$"]).assert_eq(
-        &get_latex(&vec![(0_u8.into(), Complex64::new(-1.0, 0.0))], 2)
+        &get_state_latex(&vec![(0_u8.into(), Complex64::new(-1.0, 0.0))], 2)
             .expect("expected valid latex"),
     );
     expect!([r"$|\psi\rangle = -i|00\rangle$"]).assert_eq(
-        &get_latex(&vec![(0_u8.into(), Complex64::new(0.0, -1.0))], 2)
+        &get_state_latex(&vec![(0_u8.into(), Complex64::new(0.0, -1.0))], 2)
             .expect("expected valid latex"),
     );
     expect!([r"$|\psi\rangle =  e^{-2 i \pi / 3}|00\rangle$"]).assert_eq(
-        &get_latex(
+        &get_state_latex(
             &vec![(
                 0_u8.into(),
                 Complex64::new((-2.0 * PI / 3.0).cos(), (-2.0 * PI / 3.0).sin()),
@@ -663,7 +730,7 @@ fn check_get_latex() {
         .expect("expected valid latex"),
     );
     expect!([r"$|\psi\rangle = \left( 1+\frac{\sqrt{2}}{2}i \right)|00\rangle+\left( 1+\frac{\sqrt{2}}{2}i \right)|10\rangle$"])
-    .assert_eq(&get_latex(
+    .assert_eq(&get_state_latex(
         &vec![
             (0_u8.into(), Complex64::new(1.0, 1.0 / 2.0_f64.sqrt())),
             (2_u8.into(), Complex64::new(1.0, 1.0 / 2.0_f64.sqrt())),
@@ -672,6 +739,52 @@ fn check_get_latex() {
     ).expect("expected valid latex"));
 }
 
+#[test]
+fn check_get_matrix_latex() {
+    expect!([r#"$ \begin{bmatrix} 0 & 1 \\ i & \left( 1+i \right) \\ \end{bmatrix} $"#]).assert_eq(
+        &get_matrix_latex(&vec![
+            vec![Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
+            vec![Complex64::new(0.0, 1.0), Complex64::new(1.0, 1.0)],
+        ]),
+    );
+    expect!([r#"$ \begin{bmatrix} -\left( 1-i \right) & -1 \\ -i & -\left( 1+i \right) \\ \end{bmatrix} $"#]).assert_eq(
+        &get_matrix_latex(&vec![
+            vec![Complex64::new(-1.0, 1.0), Complex64::new(-1.0, 0.0)],
+            vec![Complex64::new(0.0, -1.0), Complex64::new(-1.0, -1.0)],
+        ]),
+    );
+    expect!([r#"$ \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{i}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}} & -\frac{i}{\sqrt{2}} \\ \end{bmatrix} $"#]).assert_eq(&get_matrix_latex(&vec![
+        vec![
+            Complex64::new(FRAC_1_SQRT_2, 0.0),
+            Complex64::new(0.0, FRAC_1_SQRT_2),
+        ],
+        vec![
+            Complex64::new(-FRAC_1_SQRT_2, 0.0),
+            Complex64::new(0.0, -FRAC_1_SQRT_2),
+        ],
+    ]));
+    expect!([r#"$ \begin{bmatrix} \frac{1}{2} & \frac{i}{2} \\ -\frac{1}{2} & -\frac{i}{2} \\ \end{bmatrix} $"#]).assert_eq(&get_matrix_latex(&vec![
+        vec![
+            Complex64::new(0.5, 0.0),
+            Complex64::new(0.0, 0.5),
+        ],
+        vec![
+            Complex64::new(-0.5, 0.0),
+            Complex64::new(0.0, -0.5),
+        ],
+    ]));
+    expect!([r#"$ \begin{bmatrix} \frac{1}{2} + \frac{i}{2} & -\frac{1}{2} - \frac{i}{2} \\ -\frac{1}{2} + \frac{i}{2} & \frac{1}{2} - \frac{i}{2} \\ \end{bmatrix} $"#]).assert_eq(&get_matrix_latex(&vec![
+        vec![
+            Complex64::new(0.5, 0.5),
+            Complex64::new(-0.5, -0.5),
+        ],
+        vec![
+            Complex64::new(-0.5, 0.5),
+            Complex64::new(0.5, -0.5),
+        ],
+    ]));
+}
+
 #[test]
 fn check_get_latex_perf() {
     // This is not a CI gate for performance, just prints out data.
@@ -689,7 +802,7 @@ fn check_get_latex_perf() {
     ];
 
     expect!([r"$|\psi\rangle = \frac{1}{2}|00\rangle+\frac{1}{2} e^{ i \pi / 4}|01\rangle+\frac{1}{2}i|10\rangle+\frac{1}{2} e^{3 i \pi / 4}|11\rangle$"])
-    .assert_eq(&get_latex(
+    .assert_eq(&get_state_latex(
         &state,
         2,
     ).expect("expected valid latex"));
@@ -698,7 +811,7 @@ fn check_get_latex_perf() {
     let start = Instant::now();
     let mut l: usize = 0;
     for _ in 0..1_000 {
-        let s = get_latex(&state, 2);
+        let s = get_state_latex(&state, 2);
         l += s.map_or(0, |s| s.len());
     }
     println!(
diff --git a/compiler/qsc_eval/src/tests.rs b/compiler/qsc_eval/src/tests.rs
index e10fe65a5c..d2e183f684 100644
--- a/compiler/qsc_eval/src/tests.rs
+++ b/compiler/qsc_eval/src/tests.rs
@@ -3724,7 +3724,7 @@ fn controlled_operation_with_duplicate_controls_fails() {
                                 1,
                             ),
                             item: LocalItemId(
-                                132,
+                                134,
                             ),
                         },
                         caller: PackageId(
@@ -3774,7 +3774,7 @@ fn controlled_operation_with_target_in_controls_fails() {
                                 1,
                             ),
                             item: LocalItemId(
-                                132,
+                                134,
                             ),
                         },
                         caller: PackageId(
diff --git a/compiler/qsc_partial_eval/src/lib.rs b/compiler/qsc_partial_eval/src/lib.rs
index 3f29fffbb3..fd20e3dbf2 100644
--- a/compiler/qsc_partial_eval/src/lib.rs
+++ b/compiler/qsc_partial_eval/src/lib.rs
@@ -1304,6 +1304,7 @@ impl<'a> PartialEvaluator<'a> {
             // The following intrinsic operations and functions are no-ops.
             "BeginEstimateCaching" => Ok(Value::Bool(true)),
             "DumpRegister"
+            | "DumpOperation"
             | "AccountForEstimatesInternal"
             | "BeginRepeatEstimatesInternal"
             | "EndRepeatEstimatesInternal"
diff --git a/library/src/tests/diagnostics.rs b/library/src/tests/diagnostics.rs
index 449cd25a1e..a7789488c1 100644
--- a/library/src/tests/diagnostics.rs
+++ b/library/src/tests/diagnostics.rs
@@ -2,6 +2,7 @@
 // Licensed under the MIT License.
 
 use super::test_expression;
+use expect_test::expect;
 use qsc::interpret::Value;
 
 #[test]
@@ -41,6 +42,95 @@ fn check_operations_are_equal() {
     );
 }
 
+#[test]
+fn check_dumpoperation_for_i() {
+    let output = test_expression(
+        "Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => I(qs[0]))",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        1.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 1.0000+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_for_x() {
+    let output = test_expression(
+        "Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => X(qs[0]))",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.0000+0.0000𝑖 1.0000+0.0000𝑖
+        1.0000+0.0000𝑖 0.0000+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_for_h() {
+    let output = test_expression(
+        "Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => H(qs[0]))",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.7071+0.0000𝑖 0.7071+0.0000𝑖
+        0.7071+0.0000𝑖 −0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_for_y() {
+    let output = test_expression(
+        "Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => Y(qs[0]))",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.0000+0.0000𝑖 0.0000−1.0000𝑖
+        0.0000+1.0000𝑖 0.0000+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_for_ccnot() {
+    let output = test_expression(
+        "Microsoft.Quantum.Diagnostics.DumpOperation(3, qs => CCNOT(qs[0], qs[1], qs[2]))",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖
+        0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 0.0000+0.0000𝑖 1.0000+0.0000𝑖 0.0000+0.0000𝑖
+    "#]].assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_with_extra_qubits_allocated() {
+    let output = test_expression(
+        "{use qs = Qubit[2]; Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => H(qs[0]))}",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.7071+0.0000𝑖 0.7071+0.0000𝑖
+        0.7071+0.0000𝑖 −0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
 #[test]
 fn check_start_stop_counting_operation_called_3_times() {
     test_expression(
@@ -292,3 +382,45 @@ fn check_counting_qubits_works_with_allocation_in_operation_calls() {
         &Value::Bool(true),
     );
 }
+
+#[test]
+fn check_dumpoperation_with_extra_qubits_in_superposition() {
+    let output = test_expression(
+        "{use qs = Qubit[2]; H(qs[0]); Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => H(qs[0])); Reset(qs[0]);}",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.7071+0.0000𝑖 0.7071+0.0000𝑖
+        0.7071+0.0000𝑖 −0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_with_extra_qubits_global_phase_reflected_in_matrix() {
+    let output = test_expression(
+        "{use qs = Qubit[2]; R(PauliI, Std.Math.PI() / 2.0, qs[0]); Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => H(qs[0])); Reset(qs[0]);}",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.5000−0.5000𝑖 0.5000−0.5000𝑖
+        0.5000−0.5000𝑖 −0.5000+0.5000𝑖
+    "#]]
+    .assert_eq(&output);
+}
+
+#[test]
+fn check_dumpoperation_with_extra_qubits_relative_phase_not_reflected_in_matrix() {
+    let output = test_expression(
+        "{use qs = Qubit[2]; R1(Std.Math.PI() / 2.0, qs[0]); Microsoft.Quantum.Diagnostics.DumpOperation(1, qs => H(qs[0])); Reset(qs[0]);}",
+        &Value::unit(),
+    );
+    expect![[r#"
+        MATRIX:
+        0.7071+0.0000𝑖 0.7071+0.0000𝑖
+        0.7071+0.0000𝑖 −0.7071+0.0000𝑖
+    "#]]
+    .assert_eq(&output);
+}
diff --git a/library/std/src/Std/Diagnostics.qs b/library/std/src/Std/Diagnostics.qs
index 39a96fe572..d71e87fe62 100644
--- a/library/std/src/Std/Diagnostics.qs
+++ b/library/std/src/Std/Diagnostics.qs
@@ -64,6 +64,46 @@ function DumpRegister(register : Qubit[]) : Unit {
     body intrinsic;
 }
 
+/// # Summary
+/// Given an operation, dumps the matrix representation of the operation action on the given
+/// number of qubits.
+///
+/// # Input
+/// ## nQubits
+/// The number of qubits on which the given operation acts.
+/// ## op
+/// The operation that is to be diagnosed.
+///
+/// # Remarks
+/// When run on the sparse-state simulator, the following snippet
+/// will output the matrix
+/// $\left(\begin{matrix} 0.0 & 0.707 \\\\ 0.707 & 0.0\end{matrix}\right)$:
+///
+/// ```qsharp
+/// operation DumpH() : Unit {
+///     DumpOperation(1, qs => H(qs[0]));
+/// }
+/// ```
+/// Calling this operation has no observable effect from within Q#.
+/// Note that if `DumpOperation` is called when there are other qubits allocated,
+/// the matrix displayed may reflect any global phase that has accumulated from operations
+/// on those other qubits.
+@SimulatableIntrinsic()
+operation DumpOperation(nQubits : Int, op : Qubit[] => Unit is Adj) : Unit {
+    use (targets, extra) = (Qubit[nQubits], Qubit[nQubits]);
+    for i in 0..nQubits - 1 {
+        H(targets[i]);
+        CNOT(targets[i], extra[i]);
+    }
+    op(targets);
+    DumpMatrix(targets + extra);
+    ResetAll(targets + extra);
+}
+
+function DumpMatrix(qs : Qubit[]) : Unit {
+    body intrinsic;
+}
+
 /// Checks whether a qubit is in the |0⟩ state, returning true if it is.
 ///
 /// # Description
@@ -320,6 +360,7 @@ operation StopCountingQubits() : Int {
 export
     DumpMachine,
     DumpRegister,
+    DumpOperation,
     CheckZero,
     CheckAllZero,
     Fact,
diff --git a/npm/qsharp/src/compiler/common.ts b/npm/qsharp/src/compiler/common.ts
index 1656e5f3a5..9522c1c7ab 100644
--- a/npm/qsharp/src/compiler/common.ts
+++ b/npm/qsharp/src/compiler/common.ts
@@ -19,6 +19,12 @@ interface DumpMsg {
   stateLatex: string | null;
 }
 
+interface MatrixMsg {
+  type: "Matrix";
+  matrix: number[][][]; // Array or rows, which are an array of elements, which are complex numbers as [re, im]
+  matrixLatex: string;
+}
+
 interface MessageMsg {
   type: "Message";
   message: string;
@@ -29,7 +35,7 @@ interface ResultMsg {
   result: Result;
 }
 
-type EventMsg = ResultMsg | DumpMsg | MessageMsg;
+type EventMsg = ResultMsg | DumpMsg | MatrixMsg | MessageMsg;
 
 function outputAsResult(msg: string): ResultMsg | null {
   try {
@@ -73,12 +79,29 @@ function outputAsDump(msg: string): DumpMsg | null {
   return null;
 }
 
+function outputAsMatrix(msg: string): MatrixMsg | null {
+  try {
+    const obj = JSON.parse(msg);
+    if (obj?.type == "Matrix" && Array.isArray(obj.matrix)) {
+      return obj as MatrixMsg;
+    }
+  } catch {
+    return null;
+  }
+  return null;
+}
+
 export function eventStringToMsg(msg: string): EventMsg | null {
-  return outputAsResult(msg) || outputAsMessage(msg) || outputAsDump(msg);
+  return (
+    outputAsResult(msg) ||
+    outputAsMessage(msg) ||
+    outputAsDump(msg) ||
+    outputAsMatrix(msg)
+  );
 }
 
 export type ShotResult = {
   success: boolean;
   result: string | VSDiagnostic;
-  events: Array<MessageMsg | DumpMsg>;
+  events: Array<MessageMsg | DumpMsg | MatrixMsg>;
 };
diff --git a/npm/qsharp/src/compiler/compiler.ts b/npm/qsharp/src/compiler/compiler.ts
index a66e8c0642..39d79432b6 100644
--- a/npm/qsharp/src/compiler/compiler.ts
+++ b/npm/qsharp/src/compiler/compiler.ts
@@ -273,6 +273,12 @@ export function onCompilerEvent(msg: string, eventTarget: IQscEventTarget) {
     case "Result":
       qscEvent = makeEvent("Result", qscMsg.result);
       break;
+    case "Matrix":
+      qscEvent = makeEvent("Matrix", {
+        matrix: qscMsg.matrix,
+        matrixLatex: qscMsg.matrixLatex,
+      });
+      break;
     default:
       log.never(msgType);
       throw "Unexpected message type";
@@ -295,5 +301,5 @@ export const compilerProtocol: ServiceProtocol<ICompiler, QscEventData> = {
     run: "requestWithProgress",
     checkExerciseSolution: "requestWithProgress",
   },
-  eventNames: ["DumpMachine", "Message", "Result"],
+  eventNames: ["DumpMachine", "Matrix", "Message", "Result"],
 };
diff --git a/npm/qsharp/src/compiler/events.ts b/npm/qsharp/src/compiler/events.ts
index 20c21e23fe..7b65dec6e0 100644
--- a/npm/qsharp/src/compiler/events.ts
+++ b/npm/qsharp/src/compiler/events.ts
@@ -9,6 +9,7 @@ import { IServiceEventTarget } from "../workers/common.js";
 export type QscEventData =
   | { type: "Message"; detail: string }
   | { type: "DumpMachine"; detail: { state: Dump; stateLatex: string | null } }
+  | { type: "Matrix"; detail: { matrix: number[][][]; matrixLatex: string } }
   | { type: "Result"; detail: Result };
 
 export type QscEvents = Event & QscEventData;
@@ -81,6 +82,7 @@ export class QscEventTarget implements IQscEventTarget {
         this.onDumpMachine(ev.detail),
       );
       this.addEventListener("Result", (ev) => this.onResult(ev.detail));
+      this.addEventListener("Matrix", (ev) => this.onMatrix(ev.detail));
     }
   }
 
@@ -93,6 +95,19 @@ export class QscEventTarget implements IQscEventTarget {
     this.queueUiRefresh();
   }
 
+  private onMatrix(detail: { matrix: number[][][]; matrixLatex: string }) {
+    this.ensureActiveShot();
+
+    const shotIdx = this.results.length - 1;
+    this.results[shotIdx].events.push({
+      type: "Matrix",
+      matrix: detail.matrix,
+      matrixLatex: detail.matrixLatex,
+    });
+
+    this.queueUiRefresh();
+  }
+
   private onDumpMachine(detail: { state: Dump; stateLatex: string | null }) {
     this.ensureActiveShot();
 
diff --git a/npm/qsharp/src/debug-service/debug-service.ts b/npm/qsharp/src/debug-service/debug-service.ts
index dee78912f6..1f40daf0d5 100644
--- a/npm/qsharp/src/debug-service/debug-service.ts
+++ b/npm/qsharp/src/debug-service/debug-service.ts
@@ -163,6 +163,12 @@ export function onCompilerEvent(msg: string, eventTarget: IQscEventTarget) {
     case "Result":
       qscEvent = makeEvent("Result", qscMsg.result);
       break;
+    case "Matrix":
+      qscEvent = makeEvent("Matrix", {
+        matrix: qscMsg.matrix,
+        matrixLatex: qscMsg.matrixLatex,
+      });
+      break;
     default:
       log.never(msgType);
       throw "Unexpected message type";
@@ -190,5 +196,5 @@ export const debugServiceProtocol: ServiceProtocol<
     evalStepOut: "requestWithProgress",
     dispose: "request",
   },
-  eventNames: ["DumpMachine", "Message", "Result"],
+  eventNames: ["DumpMachine", "Message", "Matrix", "Result"],
 };
diff --git a/pip/qsharp/_native.pyi b/pip/qsharp/_native.pyi
index 97492c78e1..8f8b811996 100644
--- a/pip/qsharp/_native.pyi
+++ b/pip/qsharp/_native.pyi
@@ -103,40 +103,6 @@ class TargetProfile(Enum):
     Describes the unrestricted set of capabilities required to run any Q# program.
     """
 
-class StateDumpData:
-    """
-    A state dump returned from the Q# interpreter.
-    """
-
-    qubit_count: int
-    """
-    The number of allocated qubits at the time of the dump.
-    """
-
-    def get_dict(self) -> dict:
-        """
-        Get the amplitudes of the state vector as a dictionary from state integer to
-        complex amplitudes.
-        """
-        ...
-
-    def __repr__(self) -> str: ...
-    def __str__(self) -> str: ...
-    def _repr_html_(self) -> str: ...
-    def _repr_latex_(self) -> Optional[str]: ...
-
-class Output:
-    """
-    An output returned from the Q# interpreter.
-    Outputs can be a state dumps or messages. These are normally printed to the console.
-    """
-
-    def __repr__(self) -> str: ...
-    def __str__(self) -> str: ...
-    def _repr_html_(self) -> str: ...
-    def _repr_latex_(self) -> Optional[str]: ...
-    def state_dump(self) -> Optional[StateDumpData]: ...
-
 class Interpreter:
     """A Q# interpreter."""
 
@@ -279,6 +245,37 @@ class Pauli(Enum):
     Y: int
     Z: int
 
+class Output:
+    """
+    An output returned from the Q# interpreter.
+    Outputs can be a state dumps or messages. These are normally printed to the console.
+    """
+
+    def __repr__(self) -> str: ...
+    def __str__(self) -> str: ...
+    def _repr_markdown_(self) -> Optional[str]: ...
+    def state_dump(self) -> Optional[StateDumpData]: ...
+
+class StateDumpData:
+    """
+    A state dump returned from the Q# interpreter.
+    """
+
+    """
+    The number of allocated qubits at the time of the dump.
+    """
+    qubit_count: int
+
+    """
+    Get the amplitudes of the state vector as a dictionary from state integer to
+    complex amplitudes.
+    """
+    def get_dict(self) -> dict: ...
+    def __repr__(self) -> str: ...
+    def __str__(self) -> str: ...
+    def _repr_markdown_(self) -> str: ...
+    def _repr_latex_(self) -> Optional[str]: ...
+
 class Circuit:
     def json(self) -> str: ...
     def __repr__(self) -> str: ...
diff --git a/pip/qsharp/_qsharp.py b/pip/qsharp/_qsharp.py
index 18ca28ed50..f617384f31 100644
--- a/pip/qsharp/_qsharp.py
+++ b/pip/qsharp/_qsharp.py
@@ -25,6 +25,16 @@
 
 _interpreter = None
 
+# Check if we are running in a Jupyter notebook to use the IPython display function
+_in_jupyter = False
+try:
+    from IPython.display import display
+
+    if get_ipython().__class__.__name__ == "ZMQInteractiveShell":  # type: ignore
+        _in_jupyter = True  # Jupyter notebook or qtconsole
+except:
+    pass
+
 
 # Reporting execution time during IPython cells requires that IPython
 # gets pinged to ensure it understands the cell is active. This is done by
@@ -177,6 +187,13 @@ def eval(source: str) -> Any:
     ipython_helper()
 
     def callback(output: Output) -> None:
+        if _in_jupyter:
+            try:
+                display(output)
+                return
+            except:
+                # If IPython is not available, fall back to printing the output
+                pass
         print(output, flush=True)
 
     return get_interpreter().interpret(source, callback)
@@ -220,6 +237,13 @@ def run(
     results: List[ShotResult] = []
 
     def print_output(output: Output) -> None:
+        if _in_jupyter:
+            try:
+                display(output)
+                return
+            except:
+                # If IPython is not available, fall back to printing the output
+                pass
         print(output, flush=True)
 
     def on_save_events(output: Output) -> None:
@@ -405,11 +429,8 @@ def __repr__(self) -> str:
     def __str__(self) -> str:
         return self.__data.__str__()
 
-    def _repr_html_(self) -> str:
-        return self.__data._repr_html_()
-
-    def _repr_latex_(self) -> Optional[str]:
-        return self.__data._repr_latex_()
+    def _repr_markdown_(self) -> str:
+        return self.__data._repr_markdown_()
 
     def check_eq(
         self, state: Union[Dict[int, complex], List[complex]], tolerance: float = 1e-10
diff --git a/pip/src/displayable_output.rs b/pip/src/displayable_output.rs
index 8edadc75d5..81456fb611 100644
--- a/pip/src/displayable_output.rs
+++ b/pip/src/displayable_output.rs
@@ -6,11 +6,15 @@ mod tests;
 
 use num_bigint::BigUint;
 use num_complex::{Complex64, ComplexFloat};
-use qsc::{fmt_basis_state_label, fmt_complex, format_state_id, get_latex, get_phase};
+use qsc::{
+    fmt_basis_state_label, fmt_complex, format_state_id, get_matrix_latex, get_phase,
+    get_state_latex,
+};
 use std::fmt::Write;
 
 #[derive(Clone)]
 pub struct DisplayableState(pub Vec<(BigUint, Complex64)>, pub usize);
+pub struct DisplayableMatrix(pub Vec<Vec<Complex64>>);
 
 impl DisplayableState {
     pub fn to_plain(&self) -> String {
@@ -53,11 +57,35 @@ impl DisplayableState {
     }
 
     pub fn to_latex(&self) -> Option<String> {
-        get_latex(&self.0, self.1)
+        get_state_latex(&self.0, self.1)
+    }
+}
+
+impl DisplayableMatrix {
+    pub fn to_plain(&self) -> String {
+        format!(
+            "MATRIX:{}",
+            self.0.iter().fold(String::new(), |mut output, row| {
+                let _ = write!(
+                    output,
+                    "\n{}",
+                    row.iter().fold(String::new(), |mut row_output, element| {
+                        let _ = write!(row_output, " {}", fmt_complex(element));
+                        row_output
+                    })
+                );
+                output
+            })
+        )
+    }
+
+    pub fn to_latex(&self) -> String {
+        get_matrix_latex(&self.0)
     }
 }
 
 pub enum DisplayableOutput {
     State(DisplayableState),
     Message(String),
+    Matrix(DisplayableMatrix),
 }
diff --git a/pip/src/interpreter.rs b/pip/src/interpreter.rs
index 780897ad01..bc7c087c9e 100644
--- a/pip/src/interpreter.rs
+++ b/pip/src/interpreter.rs
@@ -2,7 +2,7 @@
 // Licensed under the MIT License.
 
 use crate::{
-    displayable_output::{DisplayableOutput, DisplayableState},
+    displayable_output::{DisplayableMatrix, DisplayableOutput, DisplayableState},
     fs::file_system,
     interop::{
         compile_qasm3_to_qir, compile_qasm3_to_qsharp, compile_qasm_enriching_errors,
@@ -547,6 +547,7 @@ impl Output {
     fn __repr__(&self) -> String {
         match &self.0 {
             DisplayableOutput::State(state) => state.to_plain(),
+            DisplayableOutput::Matrix(matrix) => matrix.to_plain(),
             DisplayableOutput::Message(msg) => msg.clone(),
         }
     }
@@ -555,24 +556,25 @@ impl Output {
         self.__repr__()
     }
 
-    fn _repr_html_(&self) -> String {
+    fn _repr_markdown_(&self) -> Option<String> {
         match &self.0 {
-            DisplayableOutput::State(state) => state.to_html(),
-            DisplayableOutput::Message(msg) => format!("<p>{msg}</p>"),
-        }
-    }
-
-    fn _repr_latex_(&self) -> Option<String> {
-        match &self.0 {
-            DisplayableOutput::State(state) => state.to_latex(),
+            DisplayableOutput::State(state) => {
+                let latex = if let Some(latex) = state.to_latex() {
+                    format!("\n\n{latex}")
+                } else {
+                    String::default()
+                };
+                Some(format!("{}{latex}", state.to_html()))
+            }
             DisplayableOutput::Message(_) => None,
+            DisplayableOutput::Matrix(matrix) => Some(matrix.to_latex()),
         }
     }
 
     fn state_dump(&self) -> Option<StateDumpData> {
         match &self.0 {
             DisplayableOutput::State(state) => Some(StateDumpData(state.clone())),
-            DisplayableOutput::Message(_) => None,
+            DisplayableOutput::Matrix(_) | DisplayableOutput::Message(_) => None,
         }
     }
 }
@@ -619,8 +621,13 @@ impl StateDumpData {
         self.__repr__()
     }
 
-    fn _repr_html_(&self) -> String {
-        self.0.to_html()
+    fn _repr_markdown_(&self) -> String {
+        let latex = if let Some(latex) = self.0.to_latex() {
+            format!("\n\n{latex}")
+        } else {
+            String::default()
+        };
+        format!("{}{latex}", self.0.to_html())
     }
 
     fn _repr_latex_(&self) -> Option<String> {
@@ -737,6 +744,22 @@ impl Receiver for OptionalCallbackReceiver<'_> {
         Ok(())
     }
 
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> std::result::Result<(), Error> {
+        if let Some(callback) = &self.callback {
+            let out = DisplayableOutput::Matrix(DisplayableMatrix(matrix));
+            callback
+                .call1(
+                    self.py,
+                    PyTuple::new_bound(
+                        self.py,
+                        &[Py::new(self.py, Output(out)).expect("should be able to create output")],
+                    ),
+                )
+                .map_err(|_| Error)?;
+        }
+        Ok(())
+    }
+
     fn message(&mut self, msg: &str) -> core::result::Result<(), Error> {
         if let Some(callback) = &self.callback {
             let out = DisplayableOutput::Message(msg.to_owned());
diff --git a/pip/src/state_header_template.html b/pip/src/state_header_template.html
index 2093cde699..77044e1c96 100644
--- a/pip/src/state_header_template.html
+++ b/pip/src/state_header_template.html
@@ -1,10 +1,38 @@
-<table>
+<table class="qs-stateTable">
+  <style>
+    .qs-stateTable thead tr {{
+      background-color: var(
+        --vscode-list-hoverBackground,
+        var(--jp-layout-color1, inherit)
+      );
+    }}
+    .qs-stateTable th {{
+      text-align: left;
+      border: none;
+    }}
+    .qs-stateTable tbody {{
+      pointer-events: none;
+    }}
+    .qs-stateTable tbody td {{
+      text-align: left;
+      border: none;
+    }}
+    .qs-stateTable tbody td span {{
+      display: inline-block;
+    }}
+    .qs-stateTable tbody tr:nth-child(even) {{
+      background-color: var(
+        --vscode-list-hoverBackground,
+        var(--jp-layout-color1, inherit)
+      );
+    }}
+  </style>
   <thead>
     <tr>
-      <th style="text-align: left">Basis State<br />(|𝜓₁…𝜓ₙ⟩)</th>
-      <th style="text-align: left">Amplitude</th>
-      <th style="text-align: left">Measurement Probability</th>
-      <th style="text-align: left" colspan="2">Phase</th>
+      <th>Basis State<br />(|𝜓₁…𝜓ₙ⟩)</th>
+      <th>Amplitude</th>
+      <th>Measurement Probability</th>
+      <th colspan="2">Phase</th>
     </tr>
   </thead>
   <tbody>
diff --git a/pip/src/state_row_template.html b/pip/src/state_row_template.html
index 41c65c245e..19a5e95626 100644
--- a/pip/src/state_row_template.html
+++ b/pip/src/state_row_template.html
@@ -1,16 +1,16 @@
 <tr>
-  <td style="text-align: left">
-    <span style="display: inline-block">|{}⟩</span>
+  <td>
+    <span>|{}⟩</span>
   </td>
-  <td style="text-align: left">
-    <span style="display: inline-block">{}</span>
+  <td>
+    <span>{}</span>
   </td>
-  <td style="text-align: left">
+  <td>
     <progress max="100" value="{}"></progress>
-    <span style="display: inline-block">{:.4}%</span>
+    <span>{:.4}%</span>
   </td>
-  <td style="text-align: left; transform: rotate({:.4}rad)">↑</td>
-  <td style="text-align: left">
-    <span style="display: inline-block">{:.4}</span>
+  <td style="transform: rotate({:.4}rad)">↑</td>
+  <td>
+    <span>{:.4}</span>
   </td>
 </tr>
diff --git a/playground/src/results.tsx b/playground/src/results.tsx
index 5fd9cfc5a2..42a9feac01 100644
--- a/playground/src/results.tsx
+++ b/playground/src/results.tsx
@@ -4,7 +4,7 @@
 import { QscEventTarget, ShotResult, VSDiagnostic } from "qsharp-lang";
 import { useEffect, useState } from "preact/hooks";
 
-import { Histogram } from "qsharp-lang/ux";
+import { Histogram, Markdown } from "qsharp-lang/ux";
 import { StateTable } from "./state.js";
 import { ActiveTab } from "./main.js";
 
@@ -240,13 +240,15 @@ export function ResultsTab(props: {
             {resultState.currResult?.events.map((evt) => {
               return evt.type === "Message" ? (
                 <div class="message-output">&gt; {evt.message}</div>
-              ) : (
+              ) : evt.type === "DumpMachine" ? (
                 <div>
                   <StateTable
                     dump={evt.state}
                     latexDump={evt.stateLatex}
                   ></StateTable>
                 </div>
+              ) : (
+                <Markdown markdown={evt.matrixLatex} />
               );
             })}
           </div>
diff --git a/vscode/src/debugger/output.ts b/vscode/src/debugger/output.ts
index 175468b8a9..0cb478f54c 100644
--- a/vscode/src/debugger/output.ts
+++ b/vscode/src/debugger/output.ts
@@ -3,23 +3,23 @@
 
 import { QscEventTarget } from "qsharp-lang";
 
+function formatComplex(real: number, imag: number) {
+  // Format -0 as 0
+  // Also using Unicode Minus Sign instead of ASCII Hyphen-Minus
+  // and Unicode Mathematical Italic Small I instead of ASCII i.
+  const r = `${real <= -0.00005 ? "−" : " "}${Math.abs(real).toFixed(4)}`;
+  const i = `${imag <= -0.00005 ? "−" : "+"}${Math.abs(imag).toFixed(4)}𝑖`;
+  return `${r}${i}`;
+}
+
 export function createDebugConsoleEventTarget(out: (message: string) => void) {
   const eventTarget = new QscEventTarget(false);
 
   eventTarget.addEventListener("Message", (evt) => {
-    out(evt.detail);
+    out(evt.detail + "\n");
   });
 
   eventTarget.addEventListener("DumpMachine", (evt) => {
-    function formatComplex(real: number, imag: number) {
-      // Format -0 as 0
-      // Also using Unicode Minus Sign instead of ASCII Hyphen-Minus
-      // and Unicode Mathematical Italic Small I instead of ASCII i.
-      const r = `${real <= -0.00005 ? "−" : ""}${Math.abs(real).toFixed(4)}`;
-      const i = `${imag <= -0.00005 ? "−" : "+"}${Math.abs(imag).toFixed(4)}𝑖`;
-      return `${r}${i}`;
-    }
-
     function formatProbabilityPercent(real: number, imag: number) {
       const probabilityPercent = (real * real + imag * imag) * 100;
       return `${probabilityPercent.toFixed(4)}%`;
@@ -38,8 +38,7 @@ export function createDebugConsoleEventTarget(out: (message: string) => void) {
     );
     const basis = "Basis".padEnd(basisColumnWidth);
 
-    let out_str = "\n";
-    out_str += "DumpMachine:\n\n";
+    let out_str = "";
     out_str += ` ${basis} | Amplitude      | Probability | Phase\n`;
     out_str +=
       " ".padEnd(basisColumnWidth, "-") +
@@ -58,8 +57,19 @@ export function createDebugConsoleEventTarget(out: (message: string) => void) {
     out(out_str);
   });
 
+  eventTarget.addEventListener("Matrix", (evt) => {
+    const out_str = evt.detail.matrix
+      .map((row) =>
+        row.map((entry) => formatComplex(entry[0], entry[1])).join(", "),
+      )
+      .join("\n");
+
+    out(out_str + "\n");
+  });
+
   eventTarget.addEventListener("Result", (evt) => {
-    out(`\n${evt.detail.value}`);
+    out(`${evt.detail.value}`);
   });
+
   return eventTarget;
 }
diff --git a/wasm/src/lib.rs b/wasm/src/lib.rs
index fa9cbc09e6..04b425b903 100644
--- a/wasm/src/lib.rs
+++ b/wasm/src/lib.rs
@@ -12,7 +12,7 @@ use num_complex::Complex64;
 use project_system::{into_qsc_args, ProgramConfig};
 use qsc::{
     compile::{self, Dependencies},
-    format_state_id, get_latex,
+    format_state_id, get_matrix_latex, get_state_latex,
     hir::PackageId,
     interpret::{
         self,
@@ -261,7 +261,7 @@ where
         )
         .expect("writing to string should succeed");
 
-        let json_latex = serde_json::to_string(&get_latex(&state, qubit_count))
+        let json_latex = serde_json::to_string(&get_state_latex(&state, qubit_count))
             .expect("serialization should succeed");
         write!(dump_json, r#" "stateLatex": {json_latex} }} "#)
             .expect("writing to string should succeed");
@@ -269,6 +269,45 @@ where
         Ok(())
     }
 
+    fn matrix(&mut self, matrix: Vec<Vec<Complex64>>) -> Result<(), output::Error> {
+        let mut dump_json = String::new();
+
+        // Write the type and open the array or rows.
+        write!(dump_json, r#"{{"type": "Matrix","matrix": ["#)
+            .expect("writing to string should succeed");
+
+        // Map each row to a string representation of the row, and join them with commas.
+        // The row is an array, and each element is a tuple formatted as "[re, im]".
+        // e.g. {"type": "Matrix", "matrix": [
+        //   [[1, 2], [3, 4], [5, 6]],
+        //   [[7, 8], [9, 10], [11, 12]]
+        // ]}
+        let row_strings = matrix
+            .iter()
+            .map(|row| {
+                let row_str = row
+                    .iter()
+                    .map(|elem| format!("[{}, {}]", elem.re, elem.im))
+                    .collect::<Vec<_>>()
+                    .join(", ");
+                format!("[{row_str}]")
+            })
+            .collect::<Vec<_>>()
+            .join(", ");
+
+        // Close the array of rows and the JSON object.
+        let latex_string = serde_json::to_string(&get_matrix_latex(&matrix))
+            .expect("serialization should succeed");
+        write!(
+            dump_json,
+            r#"{row_strings}], "matrixLatex": {latex_string} }}"#
+        )
+        .expect("writing to string should succeed");
+
+        (self.event_cb)(&dump_json);
+        Ok(())
+    }
+
     fn message(&mut self, msg: &str) -> Result<(), output::Error> {
         let msg_json = json!({"type": "Message", "message": msg});
         (self.event_cb)(&msg_json.to_string());
diff --git a/wasm/src/tests.rs b/wasm/src/tests.rs
index b7d5bd27b0..19243be6d5 100644
--- a/wasm/src/tests.rs
+++ b/wasm/src/tests.rs
@@ -129,6 +129,42 @@ fn test_message() {
     );
     assert!(result.is_ok());
 }
+
+#[test]
+fn test_matrix() {
+    let code = r"namespace Test {
+        import Microsoft.Quantum.Diagnostics.DumpOperation;
+        operation Main() : Unit {
+            DumpOperation(2, Bell);
+        }
+
+        operation Bell(q: Qubit[]) : Unit is Adj {
+            H(q[0]);
+            CX(q[0], q[1]);
+        }
+    }";
+    let expr = "Test.Main()";
+    let count = std::cell::Cell::new(0);
+    let result = run_internal(
+        SourceMap::new([("test.qs".into(), code.into())], Some(expr.into())),
+        |msg| {
+            if msg.contains("Matrix") {
+                count.set(count.get() + 1);
+                // Check the start and end of the matrix LaTeX is formatted as expected
+                assert!(msg.contains(
+                    r"$ \\begin{bmatrix} \\frac{1}{\\sqrt{2}} & 0 & \\frac{1}{\\sqrt{2}} & 0 \\\\"
+                ));
+                assert!(msg.contains(
+                    r"\\frac{1}{\\sqrt{2}} & 0 & -\\frac{1}{\\sqrt{2}} & 0 \\\\ \\end{bmatrix} $"
+                ));
+            }
+        },
+        1,
+    );
+    assert!(result.is_ok());
+    assert!(count.get() == 1);
+}
+
 #[test]
 fn message_with_escape_sequences() {
     let code = r#"namespace Sample {