Single-shot determination of quantum phases via continuous measurements

TL;DR

This paper introduces a method using continuous weak measurements to identify quantum phases in many-body systems, successfully detecting phase transitions through measurement records alone, even in highly excited states.

Contribution

It demonstrates a novel phase detection criterion based on continuous measurements, applicable to complex quantum systems, aligning well with known phase transition points.

Findings

Phase transition points are identified accurately via measurement records.

The method works in highly excited states, not just ground states.

High measurement strength induces a Zeno regime suppressing state transitions.

Abstract

We propose that weak continuous probing may be exploited to determine and define quantum phases of complex many-body systems based on the measurement record alone. We test the resulting phase criterion in numerical simulations of measurements on the Bose-Hubbard model and the quantum Ising chain. This yields a phase transition point in reasonable agreement with the quantum phase transition in the ground state of the closed system in the thermodynamic limit, despite the system being highly excited through the measurement dynamics. At high measurement strengths, the system's response enters a Zeno regime suppressing transitions between eigenstates of the measurement operator.

Tables

Table 1. Table 1: Numerical values of t in , t fin subscript 𝑡 in subscript 𝑡 fin t_{\text{in}},t_{\text{fin}} and l 𝑙 l .

Model	$t_{\text{in}}$	$t_{\text{fin}}$	$l$
Bose-Hubbard	$1.9\times {10}^5$	$2.0\times {10}^5$	$50$
Transverse-field Ising	$4.0\times {10}^4$	$5.0\times {10}^4$	$50$

Equations

$$\hat{\mathcal{H}}={\hat{\mathcal{H}}}_0+\alpha {\hat{\mathcal{H}}}_1$$

$$I(t)=2\gamma ⟨\widehat{{\mathcal{M}}_0}⟩+\sqrt{\gamma }\text{d}W/\text{d}t$$

$$\text{d}|\bar{\psi}(t)\rangle=\Bigl{[}-i\hat{\mathcal{H}}-\frac{\gamma}{2}\hat{\mathcal{M}}_{0}^{2}+I(t)\hat{\mathcal{M}}_{0}\Bigr{]}\text{d}t|\bar{\psi}(t)\rangle,$$

$$S(\omega)=(2\pi T)^{-1}\mathbb{E}\bigl{[}|\int_{0}^{T}e^{-i\omega t}I(t)\text{d}t|^{2}\bigr{]}$$

$$\begin{gathered} S(\omega )=\frac{4{\gamma }^2}{N}{\int }_{-\infty }^{\infty }\text{Tr}[{\hat{\mathcal{M}}}_0{e}^{\mathcal{L}\tau }{\hat{\mathcal{M}}}_0]e^{-i\omega \tau }\text{d}\tau \\ =\frac{8{\gamma }^2}{N}\mathrm{Re}\left[{\int }_0^{\infty }\text{Tr}[{\hat{\mathcal{M}}}_0{e}^{\mathcal{L}\tau }{\hat{\mathcal{M}}}_0]e^{-i\omega \tau }\text{d}\tau \right], \end{gathered}$$

$$\mathcal{L}[{\hat{\mathcal{M}}}_0]=-i[\hat{\mathcal{H}},{\hat{\mathcal{M}}}_0].$$

$$S_d(\omega )=\frac{8{\gamma }^2}{N}\sum _{\mathrm{Re}({\lambda }_m)<0}\frac{-\mathrm{Re}({\lambda }_m)\mathrm{Re}(t_m)+\left[\omega -\mathrm{Im}({\lambda }_m)\right]\mathrm{Im}(t_m)}{{\left[\omega -\mathrm{Im}({\lambda }_m)\right]}^2+{\left[\mathrm{Re}({\lambda }_m)\right]}^2},$$

$$S_{0}(\omega)=\frac{8\gamma^{2}}{N}\sum_{\mathrm{Re}(\lambda_{m})=0}\bigg{[}\pi\mathrm{Re}(t_{m})\delta\left(\omega-\mathrm{Im}(\lambda_{m})\right)\\ +\mathcal{P}\left(\frac{\mathrm{Im}(t_{m})}{\omega-\mathrm{Im}(\lambda_{m})}\right)\bigg{]},$$

$$\hat{\mathcal{H}}=-J\sum _{⟨j,k⟩}({\hat{b}}_j^{†}{\hat{b}}_k+{\hat{b}}_k^{†}{\hat{b}}_j)+\frac{U}{2}\sum _j{\hat{b}}_j^{†}{\hat{b}}_j({\hat{b}}_j^{†}{\hat{b}}_j-1),$$

$$f_c={\lambda }_1/N$$

$${\rho }_{ij}^{(1)}=⟨{\psi }_0|{\hat{b}}_i^{†}{\hat{b}}_j|{\psi }_0⟩,$$

$$\hat{\mathcal{H}}=-\sum _{i=1}^N{\sigma }_i^z{\sigma }_{i+1}^z-\lambda \sum _{i=1}^N{\sigma }_i^x,$$

$$\left[\hat{\mathcal{H}},{\hat{\mathcal{M}}}_0\right]=\left[\widehat{{\mathcal{H}}_0},{\hat{\mathcal{M}}}_0\right]+\alpha \left[\widehat{{\mathcal{H}}_1},{\hat{\mathcal{M}}}_0\right],$$

$$\hat{\mathcal{H}}=-\sum _{i=1}^N{\sigma }_i^z-\lambda \sum _{i=1}^N{\sigma }_i^x,$$

$$\dot{\rho }=-i[\hat{\mathcal{H}},\rho ]+\sum _i{\gamma }_i\mathcal{D}[{\hat{L}}_i]\rho (t),$$

$$\mathcal{D}[\hat{L}]\rho=\hat{L}\rho\hat{L}^{\dagger}-\frac{1}{2}\biggl{(}\hat{L}^{\dagger}\hat{L}\rho+\rho\hat{L}^{\dagger}\hat{L}\biggr{)}.$$

$$\rho (t)=e^{\mathcal{L}t}\rho (0)=\sum _i{\hat{K}}_i(t)\rho (0){\hat{K}}_i^{†}(t),$$

$$\sum _i{\hat{K}}_i^{†}(t){\hat{K}}_i(t)=1,$$

$$\begin{gathered} \mathbb{L}=-i(1\otimes \hat{\mathcal{H}}-{\hat{\mathcal{H}}}^T\otimes 1) \\ +\sum _i\frac{{\gamma }_i}{2}\left(2{\hat{L}}_i^{\ast }\otimes {\hat{L}}_i-1\otimes {\hat{L}}_i^{†}{\hat{L}}_i-{\hat{L}}_i^T{\hat{L}}_i^{\ast }\otimes 1\right), \end{gathered}$$

$$\mathbb{L}|r_m⟩⟩$$

$${\mathbb{L}}^{†}|l_m⟩⟩$$

$$\sum _m|r_m⟩⟩⟨⟨l_m|=1.$$

$$\mathcal{L}[\mathbb{1}]=-[\hat{\mathcal{H}},\mathbb{1}]+\gamma\hat{\mathcal{M}}_{0}\mathbb{1}\hat{\mathcal{M}}_{0}\\ -\frac{\gamma}{2}\biggl{(}\hat{\mathcal{M}}_{0}^{2}\mathbb{1}+\mathbb{1}\hat{\mathcal{M}}^{2}_{0}\biggr{)}=0,$$

$${\rho }^{\text{st}}=1/N$$

$${\lambda }_t[\hat{O}]={⟨\hat{O}⟩}_{\rho }\text{d}t+\frac{\text{d}{W}_t}{\sqrt{8\mathcal{K}}}.$$

$$\begin{gathered} \text{d}|\bar{\psi }\left(t\right)⟩=\{-iH\text{d}t-\mathcal{K}{\hat{O}}^2\text{d}t \\ +4\mathcal{K}\hat{O}{\lambda }_t[\hat{O}]\}|\bar{\psi }\left(t\right)⟩, \end{gathered}$$

$$I(t)\text{d}t=\frac{\gamma }{2}{\lambda }_t[\hat{O}];$$

$$\hat{O}\to 4{\hat{\mathcal{M}}}_0,\mathcal{K}\to \frac{\gamma }{32}.$$

$$\mathcal{D}[\hat{O}]\rho =\hat{O}\rho {\hat{O}}^{†}-\frac{1}{2}\left\{{\hat{O}}^{†}\hat{O},\rho \right\},$$

$$\mathcal{D}[\hat{O}]\rho =2\hat{O}\rho {\hat{O}}^{†}-\left\{{\hat{O}}^{†}\hat{O},\rho \right\}.$$

$$\text{d}\rho =-\frac{i}{\hslash }[H,\rho ]\text{d}t+2\mathcal{K}\mathcal{D}[\hat{O}]\rho \text{d}t+\underset{=\sqrt{2\mathcal{K}}\mathcal{H}[\hat{O}]\rho \text{d}{W}_t}{\underbrace{4\mathcal{K}\mathcal{H}[\hat{O}]\rho \left({\lambda }_t[\hat{O}]-{⟨\hat{O}⟩}_{\rho }\text{d}t\right)}},$$

$$\begin{array}{rl}\rho (t+\text{d}t)& =\rho (t)+\left(-\frac{i}{\hslash }[H,\rho ]+\gamma \mathcal{D}[{\hat{\mathcal{M}}}_0]\rho \right)\text{d}t+\sqrt{\gamma }\mathcal{H}[{\hat{\mathcal{M}}}_0]\rho \text{d}{W}_t\\ & \approx \rho (t)+\sqrt{\gamma }\left({\hat{\mathcal{M}}}_0\rho (t)+\rho (t){\hat{\mathcal{M}}}_0^{†}\right)\text{d}{W}_t-\sqrt{\gamma }⟨{\hat{\mathcal{M}}}_0+{\hat{\mathcal{M}}}_0^{†}⟩\text{d}{W}_t,\end{array}$$

Full text

I Introduction

Quantum phases allow descriptions of complex systems in simpler terms than a microscopic description [1]. Distinct phases span wide areas in parameter space characterized by the fundamental excitations, which govern the system’s equilibrium properties and response to perturbations. They are used to characterize a wide range of physical phenomena, including electronic, magnetic and optical properties of solid state systems [2, 3], nuclear physics [4] and cosmological topological defects [5, 6].

We consider a closed system described by a Hamiltonian

[TABLE]

that exhibits a single quantum phase transition with an abrupt change in the order parameter at the critical value $\alpha=\alpha_{c}$. The change is often due to an avoided crossing in the energy spectrum, and the location is uncovered by a change of the ground state expectation value of an appropriate order parameter that becomes infinitely sharp in the thermodynamic limit.

In dynamical phase transitions [7, 8, 9, 10, 11, 12], the situation is considerably richer than the exploration of ground state properties since the entire spectrum contributes to the dynamics. This leads to intricate questions about excited state phase transitions [13, 14] and accidental dynamical phase transitions [15]. Rather than suddenly changing the Hamiltonian, another way to quench a quantum system is to perform a measurement [16, 17, 18, 19].

When studying phase transitions experimentally, a measurement typically destroys the system. One performs the experiment several times to acquire signal or statistics. Such averaging introduces variances into otherwise well-defined parameters, e.g., particle number [20]. Although modifications of complex systems by measurement have been studied [21, 22, 23, 24], the fundamental question – if the quantum phase of a complex system can be determined from the measurement record alone – remains, to the best of our knowledge, unanswered.

In this article, we propose a criterion for the determination of quantum phases based solely on the measurement record of a single experimental run. This proposal relies on continuous dispersive measurements. It is well known that even weak and continuous measurements induce a back action through the noisy measurement record that builds up to a substantial perturbation of the system [17]. Here, we exploit this to disturb the system and simultaneously record its response. Similar to dynamical phase transitions, the entire spectrum contributes to the response. The measurement strength sets the magnitude of the disturbance. At low measurement strength, we numerically demonstrate that we can extract information about the system’s phase transition.

After introducing our criterion, we apply it to the Bose-Hubbard model and the quantum Ising chain. We show that our criterion agrees reasonably with the known phase transition in the thermodynamic limit, despite the system is not in the ground state. We also demonstrate how the measurement strength itself becomes a parameter in the open system’s phase diagram revealing (potentially controllable) properties of strongly probed systems.

II Phase Determination

Let $\hat{\mathcal{M}}_{0}$ be a Hermitian operator satisfying $[\hat{\mathcal{H}_{0}},\hat{\mathcal{M}}_{0}]=0$ and [$\hat{\mathcal{H}_{1}},\hat{\mathcal{M}}_{0}]\neq 0,$ where $\hat{\mathcal{H}}_{0}$ and $\hat{\mathcal{H}}_{1}$ are the non-commuting Hamiltonians in Eq. (1). Consider a probe that dispersively measures $\hat{\mathcal{M}}_{0}$ with strength $\gamma$. The probe yields a measurement record $I(t)$ and disturbs the system through the measurement back-action. The experimental setup is shown in Fig. 1.

For concreteness, we assume a homodyne measurement signal given by

[TABLE]

where $\langle\cdot\rangle$ denotes the expectation value and $\text{d}W$ is a Wiener increment. The state of the system conditioned on the measurement outcome evolves according to the Itô stochastic Schrödinger equation (SSE)

[TABLE]

where $\hbar=1$ and $|\bar{\psi}\rangle$ denotes a non-normalized state [18, 25, 26, 27]. The first term in Eq. (3) describes the unitary evolution. The second and third terms include the dissipation associated with the measurement. This approach is experimentally appealing, since it allows extracting phase information from a single continuous measurement.

Here we do a simulation of such an experiment. We calculate the power spectral density (PSD)

[TABLE]

by numerically integrating the SSE, see Appendices E and G. In order to sample the typical behavior away from the initial state, we discard the initial part of the quantum trajectories. We divide the considered quantum trajectory into several parts and calculate the average PSD to obtain the noise average $\mathbb{E}$.

The average dynamics, on the other hand, over different Wiener increments with $\text{d}W^{2}=\text{d}t$ is given by the Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) master equation $\dot{\rho}=\mathcal{L}[\rho]=-i[\hat{\mathcal{H}},\rho]+\gamma\mathcal{D}[\hat{\mathcal{M}}_{0}]\rho$ with $\mathcal{D}[\hat{O}]\rho=\hat{O}\rho\hat{O}^{\dagger}-\frac{1}{2}\left\{\hat{O}^{\dagger}\hat{O},\rho\right\}$ [28, 29, 30, 25, 26, 27]. Our goal is to relate the phase properties of the system to the measurement signal’s autocorrelation function $F^{(1)}_{\textrm{hom}}(t,t+\tau)=\mathbb{E}[I(t)I(t+\tau)]$ where $\mathbb{E}$ denotes a classical expectation value over the noise realizations. This correlation is given by the quantum regression theorem as $F^{(1)}_{\textrm{hom}}(t,t+\tau)=2\gamma^{2}\text{ Tr}\left[\hat{\mathcal{M}}_{0}e^{\mathcal{L}\tau}\left\{\hat{\mathcal{M}}_{0},\rho^{\text{st}}\right\}\right]$, where $\rho^{\text{st}}$ is the stationary state, such that $\mathcal{L}[\rho^{\text{st}}]=0$ [31], see also Appendix C for details. The identity is always a stationary state since $\hat{\mathcal{M}}_{0}$ is Hermitian. To verify this, one replaces $\rho^{\text{st}}=\mathbb{1}/N$ in the GKLS master equation and uses $\mathcal{D}[\hat{\mathcal{M}}_{0}]\mathbb{1}=0$. Here $N$ is the dimension of the Hilbert space. This can be understood as the measurement back-action acting as an infinite temperature heat-bath in the long-time limit [32, 33].

Considering $\rho^{\text{st}}=\mathbb{1}/N$, the stationarity of the noise process, and making use of the quantum regression theorem, we obtain

[TABLE]

where $\mathrm{Re}$ is the real part. The front factor is particular to homodyne measurements [25], see also Appendices C and D for details. Since $\mathcal{D}[\hat{\mathcal{M}}_{0}]\hat{\mathcal{M}}_{0}=0$, we have

[TABLE]

After expanding $e^{\mathcal{L}\tau}$ in Eq. (5) and utilizing Eq. (6), we conclude that the PSD is determined by the commutation relations of $\hat{\mathcal{H}}$ and $\hat{\mathcal{M}}_{0}$.

Assuming $\mathcal{L}$ is diagonalizable with eigenvalues $\lambda_{m}$ and right (left) eigenmatrices $r_{m}\;(l_{m})$, Eq. (5) is decomposed as $S(\omega)=S_{d}(\omega)+S_{0}(\omega)$ with

[TABLE]

where $t_{m}=\text{Tr}[\hat{\mathcal{M}}_{0}r_{m}]\text{Tr}[l_{m}^{\dagger}\hat{\mathcal{M}}_{0}]$, $\mathrm{Im}$ the imaginary part, $\mathcal{P}$ the Cauchy principal value, and $\delta$ the Dirac-delta function, see Appendix D. Here $S_{d}(\omega)$ ($S_{0}(\omega)$) is the part from all of the decaying (decay-free) eigenvalues of $\mathcal{L}$. We observe that the peaks in the spectra will be located at $\mathrm{Im}(\lambda_{m})$. If $\mathrm{Im}(t_{m})=0$, the eigenvalue $\lambda_{m}$ will contribute with a Lorentzian to the spectrum. The eigenmatrices with non-vanishing $\mathrm{Im}(t_{m})$ give rise to non-Lorentzian contributions in the spectrum.

As we change the parameter $\alpha$, the system (1) undergoes a phase transition at $\alpha=\alpha_{c}$. For second order quantum phase transitions, this is attributed to the level crossings at $\alpha=\alpha_{c}$ [1]. If $\gamma=0$, the unperturbed $\mathcal{L}$ has eigenvalues and vectors $\lambda_{ij}=-i(E_{i}-E_{j})=-i\omega_{ij}$ and $r_{m}=l_{m}=|\psi_{i}\rangle\langle\psi_{j}|,$ where $\hat{\mathcal{H}}|\psi_{i}\rangle=E_{i}|\psi_{i}\rangle$. For a weakly probed system, the eigenvalues and eigenvectors of $\mathcal{L}$ are obtained perturbatively. Therefore, the PSDs – which are related to the level statistics via $\lambda_{ij},r_{m},$ and $l_{m}$, cf. Eqs. (7a) and (7b) – for the two different phases are also qualitatively different. Using this, one can identify the two distinct phases in Figs. 2 and 5. Note, we have plotted the normalized PSD $\tilde{S}(\omega)=S(\omega)/\int_{-\infty}^{\infty}S(\omega)\text{d}\omega$ in the aforementioned panels for numerical convenience.

One can use two different choices of measurements – $\hat{\mathcal{M}}_{0}$ and $\hat{\mathcal{M}}_{1}$ – that commute with different parts of $\hat{\mathcal{H}}$. The PSD corresponding to one phase for the $\hat{\mathcal{M}}_{0}$ measurement is qualitatively similar to the PSD corresponding to the other phase for the $\hat{\mathcal{M}}_{1}$ measurement. In particular, we obtain from Eq. (5) that $S(\omega)\propto\delta(\omega)$ when the measured operator commutes with the Hamiltonian.

The criterion for determining a phase transition, therefore, is detecting changes in the PSD corresponding to a particular measurement $\hat{\mathcal{M}}_{i}$. Using different $\hat{\mathcal{M}}_{i}$, one can also detect multiple phase transitions. Assuming $n_{\alpha}$ phase transitions in a Hamiltonian $\hat{\mathcal{H}}(\alpha)$, one obtains representative Hamiltonians $\hat{\mathcal{H}}(\alpha_{i})$ with $i=1,2,\ldots,n_{\alpha}+1$, where $\alpha_{i}$ is a parameter value corresponding to a particular phase. We consider $n_{\alpha}+1$ distinct measuring operators satisfying $[\hat{\mathcal{H}}(\alpha_{i}),\hat{\mathcal{M}}_{i}]=0$. Note that since $[\hat{\mathcal{M}}_{i},\hat{\mathcal{H}}(\alpha)]\neq 0$, our continuous measurement scheme is not a quantum non-demolition measurement [34, 35].

In the following, we implement this scheme to study the phase transitions in an ergodic (Bose-Hubbard) system and in an integrable (transverse-field Ising chain) Hamiltonian.

III Probed Bose-Hubbard Model

The 1D Bose-Hubbard model provided the first demonstration of a quantum phase transition in ultracold atoms [36], and it is a powerful tool for the experimental study of quantum phases [37, 38, 39, 40, 41], including studies of driven-dissipative quantum systems [23, 22, 42]. The Hamiltonian reads

[TABLE]

where the bosonic field operators are expanded in Wannier functions $\hat{\Psi}(x,t)=\sum_{j}\hat{b}_{j}(t)w_{j}(x)$, and $J$ and $U$ are the hopping and the on-site interaction, respectively. For $\alpha=U/J$ below the critical value, the system’s ground state exhibits long range phase-coherence and it is a superfluid. Above that critical value, the ground state features Fock-states on each site and the system is in the Mott-insulator phase.

Let us now dispersively probe this system with an optical cavity field aligned with the trapping lattice. The probe light is described as $\hat{a}(t)f_{a}({\bf x},\omega_{\mathrm{L}})e^{-i\omega_{\mathrm{L}}t}$ with $\omega_{\mathrm{L}}$ the probe frequency and $f_{a}({\bf x},\omega_{\mathrm{L}})$ the spatial mode function. Here we treat the system in 1D. For a Fabry-Pérot cavity, we have $f_{a}(x,\omega_{\mathrm{L}})\propto\cos(k_{\mathrm{L}}x)$ with $k_{\text{L}}$ being the wavenumber for the probe light.

We focus on two relevant cases, namely where the probe has twice the period of the trapping potential and when the probe and the lattice have the same periodicity, but a $\pi/2$ phase shift. In the former case, this leads to a measurement operator $\hat{\mathcal{M}}_{\mathrm{pop}}=m_{\mathrm{pop}}\sum_{j}\hat{b}_{2j}^{\dagger}\hat{b}_{2j}$, where $m_{\mathrm{pop}}$ is a constant calculated from the Wannier functions, see Appendix F. This operator commutes with the interaction term in (8) but not with the hopping. In the second case, we measure the sum over coherences, $\hat{\mathcal{M}}_{\mathrm{coh}}=m_{\mathrm{coh}}\sum_{j}\hat{b}_{j}^{\dagger}\hat{b}_{j+1}+h.c.$ (see Appendix F), which commutes with the hopping term but not the interaction.

We numerically calculate the PSDs for both $\hat{\mathcal{M}}_{\mathrm{pop}}$ and $\hat{\mathcal{M}}_{\mathrm{coh}}$. To perform the numerical integrations of the SSE (3) in Fig. 2(a,c), we considered a system with six sites and six particles. We used a system with a smaller Hilbert space – four sites and four particles – to obtain the PSDs using Eq. (5) in Fig. 2(b,d). For all the PSDs, the spectral range is rescaled to 20 in dimensionless units. The PSDs for a particular measurement – e.g., Fig. 2(a,b) – obtained from Eqs. (3) and (5) appear similar. This is because of the ergodicity of the Bose-Hubbard model.

We observe large values of the PSD at $\omega=0$ ($\omega\neq 0$) when the measurement operator and Hamiltonian is (is not) compatible with the quantum phase. The $\hat{\mathcal{M}}_{\text{coh}}$ PSDs in the superfluid (Mott-insulator) part is qualitatively similar to the Mott-insulator (superfluid) part of the $\hat{\mathcal{M}}_{\text{pop}}$ PSDs. The measurement for both operators gives the transition within the same order of magnitude, which is also in agreement with its value in the thermodynamic limit [43, 44, 45]. Additionally, we point out that the phase transition point from the PSDs are consistent with the behaviour of the order parameter (condensate fraction [46])

[TABLE]

in Fig. 3(a), where $\lambda_{1}$ is the largest eigenvalue of the single-particle density matrix $\rho^{(1)}$ and $N$ is the number of particles. The matrix elements of $\rho^{(1)}$ are given by

[TABLE]

where $|\psi_{0}\rangle$ is the ground state of the Hamiltonian (8). While comparing with the above critical value of the parameter, one needs to, however, keep in mind the significant finite size effects [46], see also Fig. 3(a). The numerical results for the thermodynamic limit [44, 45] are obtained by calculating the energy gap between the ground state and first excited state for different system sizes, and extrapolating to the infinite system.

III.1 Strong Measurement

The measured system has phase transitions defined and/or controlled by the measurement itself. That transition due to strong measurement is also witnessed by the record. The measurement strength is considered a free parameter and an additional dimension of the phase diagram, which then depends on the operator being measured. Figure 4 shows the PSD for measuring $\hat{\mathcal{M}}_{\text{coh}}$ with $\gamma\gg 1$ in the Bose-Hubbard model using Eq. (5). The figure shows that the measurement forces the system to evolve into eigenstates of the probed operator over a broad range of $U/J$. Since $\hat{\mathcal{M}}_{\text{coh}}$ commutes with the Bose-Hubbard hopping term, the PSD implies a superfluid phase throughout. This has been identified previously as a dynamical phase transition into a Zeno regime [47, 48, 49, 50]. In the given example, we demonstrate how a strong measurement of coherence turns a Mott-insulator into a superfluid. Performing strong measurements with other operators yield similar results.

IV Probed Transverse-Field Ising Chain

We now show that the change in the PSD reveals the phase transition in the transverse-field Ising chain, which is exactly solvable using the Jordan-Wigner transformation and is a paradigm for quantum phase transitions [1]. This model was implemented with trapped ions [51, 52], where a dynamical phase transition was observed [53]. The Hamiltonian is

[TABLE]

where $\sigma_{i}^{x,z}$ are Pauli operators, we use periodic boundary conditions, and $\lambda$ is a dimensionless parameter. As $\lambda$ is varied, the system exhibits a quantum phase transition at $\lambda_{c}=1$ in the thermodynamic limit from a ferromagnetic $\lambda<\lambda_{c}$ to a paramagnetic $\lambda>\lambda_{c}$ phase. We consider a homodyne measurement of the coupling $\hat{\mathcal{M}}_{\text{ZZ}}=\sum_{i=1}^{N}\sigma_{i}^{z}\sigma_{i+1}^{z}$ and transverse-field $\hat{\mathcal{M}}_{\text{X}}=\sum_{i=1}^{N}\sigma_{i}^{x}$.

The PSDs obtained by numerically integrating the SSE (3) for a system with $N=10$ and measurement operators $\hat{\mathcal{M}}_{\text{X}}$ and $\hat{\mathcal{M}}_{\text{ZZ}}$ are shown in Fig. 5(a,c), respectively. We also obtain the PSDs using Eq. (5) for a system with $N=6$ spins for the same measurement operators in Fig. 5(b,d). In order to compare the PSDs for different values of $\lambda$, we always rescale the Hamiltonian such that its spectrum spans the same frequency range (20 in dimensionless units).

Similar to the Bose-Hubbard PSDs, the qualitative nature of the PSDs change when we go from the ferromagnetic to the paramagnetic phase. We note that the $\hat{\mathcal{M}}_{\text{X}}$ PSDs in the ferromagnetic part is qualitatively similar to the paramagnetic part of the $\hat{\mathcal{M}}_{\text{ZZ}}$ PSDs. This is true for all the PSDs. The ferromagnetic part of $\hat{\mathcal{M}}_{\text{ZZ}}$ PSD obtained using Eq. (5) in Fig. 5(d) is similar to the paramagnetic $\hat{\mathcal{M}}_{\text{X}}$ PSDs.

To obtain the PSDs using the SSE (3), we start with the ground state of the Hamiltonian (1) at $t=0$. Moreover, if $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{i}]=0$ – e.g., when $\lambda=0$ $(1/\lambda=0)$ in the $\hat{\mathcal{M}}_{\text{ZZ}}$ $(\hat{\mathcal{M}}_{\text{X}})$ measurement in the transverse-field Ising model – the measurement process does not change the initial wavefunction. This leads to a flat PSD with no features. This is unlike the GKLS PSDs assuming $\rho^{\text{st}}=\mathbb{1}/N$, where $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{i}]=0$ results in $S(\omega)\propto\delta(\omega)$.

We have $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{i}]\neq 0$ for the parameter ranges considered in Figs. 2 and 5. Therefore, the measurement process is equivalent to an exploration of the phase space even if we start with an eigenstate of $\hat{\mathcal{H}}$. However, we believe that the integrability of the transverse-field Ising model is responsible for the absence of any peaks in the ferromagnetic $\hat{\mathcal{M}}_{\text{ZZ}}$ PSD Fig. 5(c) obtained using Eq. (3). Since we start with a mixed state $\propto\mathbb{1}$ while using Eq. (5), the $\hat{\mathcal{M}}_{\text{ZZ}}$ PSD still has a peak even in the ferromagnetic phase.

V Change in PSD due to the Commutation Relation $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{0}]$

In the foregoing analysis, we considered a Hamiltonian that depended on a single parameter $\alpha$. The measurement operator $\hat{\mathcal{M}}_{0}$ is chosen such that it commutes with one part of the Hamiltonian $\hat{\mathcal{H}_{0}}$ while not commuting with the other: $\hat{\mathcal{H}_{1}}$. Writing the commutation relation between the Hamiltonian and the measurement operator as

[TABLE]

we observe that $\hat{\mathcal{M}}_{0}$ evolves from commuting with $\hat{\mathcal{H}}$ to not commuting as we change $\alpha$. However, it is important to note that $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{0}]\neq 0$ for the range of $\alpha$ considered in Sec. III with the identification $\alpha\equiv U/J$ and in Sec. IV with $\alpha\equiv\lambda$.

Nevertheless, one needs to be careful while discerning the changes in PSDs due to phase transitions since in the finite sized systems one expects to see some changes in the PSDs simply because of the commutation properties (e.g., $[\hat{\mathcal{H}},\hat{\mathcal{M}}_{0}]$ being equal or unequal to zero). To illustrate this, we consider the Hamiltonian

[TABLE]

which does not go through a phase transition. The ground state of $\hat{\mathcal{H}_{0}}=-\sum_{i=1}^{N}\sigma_{i}^{z}$ is connected to the ground state of $\hat{\mathcal{H}_{1}}=-\sum_{i=1}^{N}\sigma_{i}^{x}$ by continuous rotations. This is demonstrated by the ground state expectation values of $\sum_{i}\sigma_{i}^{z}/N$ and $\sum_{i}\sigma_{i}^{x}/N$ in Figs. 6(a) and 6(b), both of which are similar to the inset of Fig. 3(b).

Similar to Sec. IV, we show the PSDs obtained for the measurements $\hat{\mathcal{M}}_{\text{Z}}=\sum_{i=1}^{N}\sigma_{i}^{z}$ and $\hat{\mathcal{M}}_{\text{X}}=\sum_{i=1}^{N}\sigma_{i}^{x}$ in Fig. 7 for the Hamiltonian (13) by using Eq. (5). Here we see a continuous change in the PSDs between $\lambda\approx 0.1$ and $\lambda\approx 1$. Interestingly, this change appears at the same interval in $\lambda$ for six as well as for four spin PSDs. The position of the peaks in $\omega$ are slightly different in the four and the six spin PSDs, whereas the heights remain almost unchanged.

We note that the changes in the PSDs due to the commutation relations are not as abrupt as the ones caused by the change in the Hamiltonian spectrum due to a phase transition. We believe that the commutation relations change the PSDs trivially compared to the changes occurring due to a phase transition, and these two types of changes in the PSDs can indeed be differentiated. However, to confirm this hypothesis peremptorily, one either needs to independently verify with an experiment or to perform numerics on a thermodynamically large system.

VI Summary and Outlook

We show that it is possible to detect phase transitions in the 1D Bose Hubbard model and the transverse-field Ising model by discerning the qualitative changes in the measurement signals of weak continuous measurements. To observe these changes, one need not prepare the state in a particular way or be confined to the ground state. We believe that this method of detecting the phase transition can be applied to various strongly interacting systems for a range of experimentally realizable measurement operators.

We have focused on the situation where the system Hamiltonian is known. In other situations of interest, this might not be the case. It will be interesting to investigate what can be deduced about a system’s Hamiltonian from measurement records. Furthermore, our criterion may be generalized to topological [54] and dynamical phase transitions [53], which have also been implemented successfully. Further exploration and, in particular, experiments will be needed to assess the broader applicability of continuous measurements as a probe of phase transitions.

Acknowledgements.

This work was supported by the Villum Foundation, the Independent Research Fund Denmark under Grant Number 8049-00074B, the Carlsberg Foundation, and ERC, H2020 grant 639560 (MECTRL). LFB would like to thank the Max Planck Institute for the Physics of Complex Systems for hospitality during visits to the institute.

Organization of the Appendices

We start by reviewing a few important properties of the Liouvillian and vectorization (the latter is also known as the Choi-Jamiłkowski isomorphism [55, 56]), which is used extensively in the following calculations. We then write the stochastic master equation keeping terms up to order $\sqrt{\text{d}t}$. In the process, we compare the notations of Refs. [25], [49], [27] and [26]. We go over the derivation for the autocorrelation function of the measurement record $F^{(1)}_{\textrm{hom}}(t,t+\tau)$. Starting from the expression for $F^{(1)}_{\textrm{hom}}$ and making use of the quantum regression theorem, we derive the expression for the PSD – Eq. (5) in the main text. We further simplify this using the Choi-Jamiłkowski isomorphism and obtain Eqs. (7a) and (7b) of the main text. We describe the numerical procedure for obtaining the PSDs in Figs. 2(b,d), 5(b,d) and 4 in the main text. We provide an expression and examples of pictorial representations of the matrix elements $M_{jk}$ necessary for constituting the measurement operator $\hat{\mathcal{M}}_{0}$. Finally, we describe the numerical integration procedure for obtaining the PSDs in Figs. 2(a,c) and 5(a,c).

Appendix A Liouvillian and Vectorization

Markovian dynamics of a linear and completely positive open quantum system can be described by a Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) master equation [28, 29, 30, 25, 26, 27]

[TABLE]

where

[TABLE]

The GKLS master equation is linear in $\rho$, which allows us to associate it with the so-called Liouvillian superoperator $\mathcal{L}$ satisfying $\partial_{t}\rho=\mathcal{L}\rho$. The superoperator $\mathcal{L}$ is trace preserving and generates the following completely positive trace preserving map $e^{\mathcal{L}t}$ describing the time evolution of the system:

[TABLE]

such that

[TABLE]

where the set of operators $\{\hat{K}_{i}\}$ are called Kraus operators. The above way of representing the completely positive trace preserving map is called the operator-sum representation.

Superoperators such as $\mathcal{L}$ act on the Liouville space $\mathcal{B}(\mathcal{H})$ consisting of all the linear operators acting on the Hilbert space. This space can itself be treated as a Hilbert space with the Hilbert-Schmidt inner-product $\langle\langle\hat{A}|\hat{B}\rangle\rangle=\text{Tr}(\hat{A}^{\dagger}\hat{B})$. We use the notation $|\rho\rangle\rangle$ for a vectorized state that is created by stacking the columns of $\rho$. In order to ease the calculations, we apply this vectorized notation here [57]. The vectorized representation of $\mathcal{L}$ is

[TABLE]

where $A^{T}$ denotes the transpose. Note that generally $\mathbb{L}$ is a non-Hermitian matrix.

In this paper, we are only concerned with diagonalizable Liouvillians. For non-diagonalizable Liouvillians, one needs to consider the Jordan normal form. Unlike the Hamiltonian, the Liouvillian is generally not Hermitian, i.e., the adjoint superoperator $\mathcal{L}^{\dagger}$ is not equal to $\mathcal{L}$. For this reason, the eigenvalues of $\mathcal{L}$ are generally complex and it has different right and left eigenstates satisfying

[TABLE]

We fix the normalization such that the left and right eigenstates are orthonormal $\langle\langle r_{m}|l_{n}\rangle\rangle=\delta_{mn}$, which is called the biorthogonality. Enumerating the eigenstates according to Eqs. (19a) and (19b), we obtain the following completeness relation:

[TABLE]

We assume that the open system dynamics are due to a continuous weak measurement of a single Hermitian operator $\hat{L}=\hat{\mathcal{M}}_{0}$. At long times, the unmonitored system will reach a steady state of $\mathcal{L}$ defined by $\mathcal{L}[\rho^{\text{st}}]=0$, i.e., a member of the kernel of the operator $\mathcal{L}$. For a Hermitian operator, we make the simple observation that

[TABLE]

which shows that $\mathbb{1}/N$ is always a stationary state where N is the dimension of the Hilbert space. In general, there is no guarantee that this is the only stationary state, but we assume

[TABLE]

for simplicity.

Appendix B Stochastic Master Equation

We write the stochastic master equation (SME) corresponding to the stochastic Schrödinger equation (SSE) considered in the main text. In the SME, we only keep the terms upto order $\sqrt{dt}$. In the process, we reconcile the derivations and notations of Refs. [25], [26], [27] and [49].

The definition of homodyne current in Refs. [49], [27] and [26] is as follows:

[TABLE]

The corresponding SSE is

[TABLE]

where $\left|\bar{\psi}\left(t\right)\right\rangle$ symbolizes the non-normalized wavefunction. We obtain the homodyne current and the SSE of the main text from Eqs. (23) and (24) as follows:

[TABLE]

While applying the prescription (25b) in Eq. (24), one does not change the wavefunction $\left|\bar{\psi}\left(t\right)\right\rangle$.

We now apply the prescription (25b) to the SME of Ref. [49]. Before we do so, we clarify the different definitions of the Lindblad superoperator appearing in different references. We list all the definitions below as

[TABLE]

​​For the rest of the discussion, we will be using the definition (26a). The SME conditioned on the random measurement outcome (23) is

[TABLE]

where $\mathcal{H}[\hat{O}]\rho=\hat{O}\rho+\rho\hat{O}^{\dagger}-\left\langle\hat{O}+\hat{O}^{\dagger}\right\rangle_{\rho}\rho$. All the references agree on the definition of $\mathcal{H}[\hat{O}]\rho$. Here, we consider the detector to be $100\%$ efficient. Also, since we are using the definition (26a), the coefficient of the second term is $2\mathcal{K}$ instead of $\mathcal{K}$ (cf. Eq. (5) of [49]).

Using the prescription (25b) we replace $\sqrt{2\mathcal{K}}$ by $\sqrt{\gamma}/4$ and $\hat{O}$ by $4\hat{\mathcal{M}}_{0}$. Additionally, we write $\left\langle\cdots\right\rangle_{\rho}$ as $\left\langle\cdots\right\rangle$ for notational convenience. Finally, we obtain

[TABLE]

where we have retained terms only upto order $\sqrt{\text{d}t}$. Later in this manuscript, we consider $\rho(t)$ is a priori known to be $\rho^{\text{st}}=\mathbb{1}/N$. In the above equation, $\rho(t+\text{d}t)$ is conditioned on the homodyne current until time $t$.

Appendix C Output Field Correlation Function

We revisit the derivation for the autocorrelation function of the measurement record $F^{(1)}_{\textrm{hom}}(t,t+\tau)=\mathbb{E}[I(t+\tau)I(t)]$. Note that we write the autocorrelation function as $F^{(1)}$. This is because of its relation to Glauber’s first-order coherence function. Here we follow Ref. [25] closely. The steps are as follows:

[TABLE]

While obtaining $\left\langle\hat{\mathcal{M}}_{0}\right\rangle(t+\tau)$, the trace is calculated with the density operator of Eq. (28) and by identifying $\text{d}t$ with $\tau$. The factorization in the last term of the last line is justified because $\rho(t)$ is given. Using similar argument we have

[TABLE]

We explicitly calculate the first term of the last line in Eq. (29) as follows:

[TABLE]

In the first line, $e^{\mathcal{L}\tau}$ provides the noise averaged time evolution between $t+\text{d}t$ and $t+\tau$. In the final line, we used the Itô rule. While expanding the superoperator $\mathcal{H}[\hat{\mathcal{M}}_{0}]$, we also used the fact that $\hat{\mathcal{M}}_{0}$ is self-adjoint. Substituting Eq. (31) into Eq. (29), we obtain

[TABLE]

To simplify the last term of Eq. (32), we note the following:

[TABLE]

where in the second to last line while performing the $\mathbb{E}$ operation, we recall that $\text{d}W_{t}$ and $\mathcal{H}[\hat{\mathcal{M}}_{0}]\rho(t)$ are statistically independent. Using Eq. (33) into Eq. (32) and substituting $\rho(t)=\rho^{\text{st}}=\mathbb{1}/N$, we finally obtain the autocorrelation function as

[TABLE]

where we have used vectorization and that $\hat{\mathcal{M}}_{0}$ is self-adjoint to write the final expression. The $\delta$-function in this formula arises due to the local oscillator shot noise or vacuum noise.

Appendix D Homodyne Spectrum: Exact Analytical Result

The power spectral density is the Fourier transformation of $F^{(1)}_{\textrm{hom}}(t,t+\tau)$ with the $\delta$-function dropped. First, note that the result of $e^{\mathcal{L}\tau}$ acting on the Hermitian operator $\hat{\mathcal{M}}_{0}\rho^{\text{st}}+\rho^{\text{st}}\hat{\mathcal{M}}^{\dagger}_{0}$ can be written using the operator-sum representation, see Eq. (16). As a result, the autocorrelation function (without the $\delta$-function) is the trace of a Hermitian operator, which is real. Moreover, the autocorrelation function is an even function in $\tau$. Using these properties, we obtain

[TABLE]

where

[TABLE]

Here we have inserted unity (20) and used that $|r_{m}\rangle\rangle$ is a right eigenstate.

For calculating the integral we must consider two special cases

[TABLE]

where $\mathcal{P}$ denotes the Cauchy principal value. We do not need to consider the case for $\mathrm{Re}(\lambda_{m})>0$ since $\mathbb{L}$ cannot have eigenvalues with positive real part. Notice the Cauchy principal value is similar to the result in Eq. (37a) but there is an additional $\delta$-function contribution in the second integral. Plugging this into Eq. (36) gives the result

[TABLE]

where $t_{m}=\text{Tr}[A^{\dagger}r_{m}]\text{Tr}[l_{m}^{\dagger}B]$. This is also the result given in the main text in Eqs. (7a) and (7b).

Appendix E Homodyne Spectrum: Numerical Computation

In this section, we describe how we compute $S(\omega)$ numerically. In principle, we could diagonalize $\mathbb{L}$ and utilize Eqs. (7a) and (7b) of the main text. However, this is not feasible due to the large dimensionality of the Liouville space. Instead, from Eq. (36) we observe that

[TABLE]

For numerical convenience, we have introduced $|\hat{\xi}\rangle\rangle$ as the solution to the linear equation system

[TABLE]

This equation must be solved for each value of $\omega$. The matrix $(i\omega\mathbb{1}-\mathbb{L})$ preserves the sparsity of the original Hamiltonian. Even then, we could only compute the numerical spectrum for the 1D Bose-Hubbard model with four sites and four particles, and for the transverse-field Ising model with $N=6$ spins. Recall that we could simulate the SSE for the 1D Bose-Hubbard model with six sites and six particles, and for the transverse-field with $N=10$ spins.

We observe from Eqs. (35) and (38) that $S(\omega)$ is singular if $\mathrm{Re}(\lambda_{m})=0$ and $\mathrm{Im}(\lambda_{m})=\omega$. The system cannot be solved for $\omega=0$, since $\mathbb{L}$ is singular. In the transverse-field Ising (Bose-Hubbard) model, we solve Eq. (40) for 204 (818) linearly spaced values of $\omega$ between 0.04 (0.01) and 8.00 in dimensionless units. We do not encounter any singularities for these frequency grids. However, in our numerical experience, the system becomes much harder to solve as $\omega\rightarrow 0$.

Appendix F The Measurement Operators in the Bose-Hubbard Model

We write a dispersive measurement operator as

[TABLE]

where $\hat{b}_{i}^{\dagger}$ creates a boson at the $i$th optical lattice site. One obtains the matrix elements in terms of the Wannier functions as

[TABLE]

where $g$ denotes the coupling strength between the probe laser and the ultracold atomic system, $f_{a}(x,\omega_{L})$ is the spatial mode function, and $\Delta$ is the detuning of the probe from the atomic transition [58, 59, 60]. For the two measurement operators considered in the main text, the entries of the matrices are displayed as images in Fig. 8. From there we have

[TABLE]

where $\textrm{mod}()$ denotes the modulo operation. We ignore the term $d_{\mathrm{coh}}\delta_{j,k}$ in our numerical integration. This is because this term leads to a constant shift $\hat{\mathcal{C}}=d_{\mathrm{coh}}N_{b}$ in $\hat{\mathcal{M}}_{\mathrm{coh}}$ with $N_{b}=6$ being the total number of bosons in the system, and the normalized Itô SSE

[TABLE]

remains unchanged under the transformation $\hat{\mathcal{M}}_{0}\rightarrow\hat{\mathcal{M}}_{0}+\hat{\mathcal{C}}$, where $\hat{\mathcal{C}}$ is a constant operator.

Appendix G Numerical Integration of the Stochastic Schrödinger Equation

We start with the Itô SSE

[TABLE]

that describes the time evolution of a non-normalized wavefunction $|\bar{\psi}\rangle$. We have written the homodyne measurement signal $I(t)$ as

[TABLE]

To obtain the PSDs in Figs. 1(b,d,f,h) of the main text, we use the Stratonovich form of the SSE [25, 61]

[TABLE]

We need this form because the chain rule for Stratonovich equations is equivalent to the chain rule of conventional calculus.

Discretizing the full time interval $(0,t_{\text{fin}}]$, we write the wavefunction at the $(j+1)$th step as

[TABLE]

where

[TABLE]

$\delta t=t_{j+1}-t_{j}$ is the infinitesimal time increment, and $S_{j}$ is a random number drawn from a standard normal distribution. To compute $e^{\hat{\mathbb{G}}(t_{j})}\left|\psi(t_{j})\right\rangle$, we use a Krylov subspace projection technique. Instead of computing the matrix exponential in isolation, this technique directly computes the action of the exponential operator on the wavefunction. Although we used the normalized wavefunction $\left|\psi(t_{j})\right\rangle$ on the right-hand side of Eq. (48), we need to normalize the wavefunction again at the $(j+1)$th step using $\left|\psi(t_{j+1})\right\rangle=|\bar{\psi}(t_{j+1})\rangle/\sqrt{\langle\bar{\psi}(t_{j+1})|\bar{\psi}(t_{j+1})\rangle}$.

After numerically obtaining the trajectories $\{\left|\psi(t_{j})\right\rangle,I(t_{j})\}$ for all the time steps in the interval $(0,t_{\text{fin}}]$, we discard the initial transients corresponding to the part $(0,t_{\text{in}}]$. To obtain a noise averaged smoother PSD, we divide the considered quantum trajectory into $l$ parts and calculate the average PSD. The values of $t_{\text{in}},\;t_{\text{fin}}$ and $l$ for Bose-Hubbard and the transverse-field Ising model are given in Table 1. We have considered the infinitesimal time increment $\delta t$ to be $0.01$ for all the trajectories. We also show how the PSD depend on $l$ in Fig. 9.