Time-dependent quantum correlations in phase space

TL;DR

This paper introduces a regularized, smooth quasi-probability distribution in phase space to visualize and analyze time-dependent quantum correlations of light, applicable even for highly singular states, and demonstrates its use in characterizing optical parametric processes.

Contribution

It develops a regularized, multi-time P functional for quantum correlations, enabling experimental access and analysis of complex quantum states in phase space.

Findings

Successfully characterizes optical parametric processes with frequency mismatch

Provides a measurement scheme for direct experimental access to the P functional

Extends the applicability of phase-space methods to highly singular quantum states

Abstract

General quasi-probabilities are introduced to visualize time-dependent quantum correlations of light in phase space. They are based on the generalization of the Glauber-Sudarshan P function to a time-dependent P functional [W. Vogel, Phys. Rev. Lett. 100, 013605 (2008)], which fully describes temporal correlations of radiation fields on the basis of continuous phase-space distributions. This approach is nontrivial as the P functional itself is highly singular for many quantum states and nonlinear processes. In general, it neither yields a well-behaved nor an experimentally accessible description of quantum stochastic processes. Our regularized version of this multi-time-dependent quasi-probability is a smooth function and applies to stronger divergences compared to the single-time and multi-mode scenario. The technique is used to characterize an optical parametric process with frequency mismatch. A measurement scheme, together with a sampling approach, is provided which yields direct experimental access to the regularized P functional from measured data.

Equations

$$\hat{\rho }=\int d^2\alpha P(\alpha )|\alpha ⟩⟨\alpha |,$$

$$P_{\Omega }(\alpha )={\mathcal{F}}_{\beta }[{\Omega }_w(\beta )\Phi (\beta )](\alpha )=(P\ast {\tilde{\Omega }}_w)(\alpha ),$$

$$P_{\Omega }(\mathbf{\alpha })={\mathcal{F}}_{\mathbf{\beta }}[{\Omega }_w(\mathbf{\beta })\Phi (\mathbf{\beta })](\mathbf{\alpha }),$$

$$|\Phi (\mathbf{\beta })|\le \exp [|\mathbf{\beta }{|}^2/2],$$

$$|\Phi (\mathbf{\beta })|=\left|⟨:\prod _{j=1}^N\hat{D}({\beta }_j):⟩\right|=\left|⟨e^{{\sum }_{j=1}^N{\beta }_j{\hat{a}}_j^{†}}{e}^{-{\sum }_{j=1}^N{\beta }_j^{\ast }{\hat{a}}_j}⟩\right|,$$

$$|\Phi (\mathbf{\beta })|=\left|⟨\prod _{j=1}^N\hat{D}({\beta }_j)⟩\right|e^{|\mathbf{\beta }{|}^2/2}.$$

$$\displaystyle\big{\langle}\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\hat{O}[\{\hat{a}(t_{i})\}_{i=1}^{k}]\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\big{\rangle}=$$

$$P[\{\alpha_{i};t_{i}\}_{i=1}^{k}]=\Big{\langle}\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\prod_{i=1}^{k}\hat{\delta}(\hat{a}(t_{i})-\alpha_{i})\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\Big{\rangle},$$

$$P[\{{\alpha }_i;t_i{\}}_{i=1}^k]={\mathcal{F}}_{\{{\beta }_i{\}}_{i=1}^k}[\Phi (\{{\beta }_i;t_i{\}}_{i=1}^k)](\{{\alpha }_i{\}}_{i=1}^k),$$

$$\Phi (\{{\beta }_i;t_i{\}}_{i=1}^k)=⟨\begin{array}{c}\circ \\ \circ \end{array}\prod _{i=1}^k\hat{D}({\beta }_i;t_i)\begin{array}{c}\circ \\ \circ \end{array}⟩$$

$${\Phi }_{\Omega }(\{{\beta }_i;t_i{\}}_{i=1}^k)=\Phi (\{{\beta }_i;t_i{\}}_{i=1}^k){\Omega }_{\mathbf{w}}(\{{\beta }_i;t_i{\}}_{i=1}^k).$$

$$P_{\Omega }[\{{\alpha }_i;t_i{\}}_{i=1}^k]={\mathcal{F}}_{\{{\beta }_i{\}}_{i=1}^k}[{\Phi }_{\Omega }(\{{\beta }_i;t_i{\}}_{i=1}^k)](\{{\alpha }_i{\}}_{i=1}^k),$$

$$\mathrm{tri}(z)=\left\{\begin{array}[]{ll}(1+z)&\text{ for }z\in[-1,0],\\ (1-z)&\text{ for }z\in[0,1],\\ 0&\text{ else. }\end{array}\right.$$

$${\Omega }_{\mathbf{w}}(\{{\beta }_i{\}}_{i=1}^k)=\prod _{i=1}^k\left(\mathrm{tri}(\mathrm{Re}[{\beta }_i]/w_i)\mathrm{tri}(\mathrm{Im}[{\beta }_i]/w_i)\right),$$

$${\hat{H}}_{\mathrm{int}}(t)=\hslash \kappa ({\hat{a}}^{†2}{e}^{-i\delta t}+{\hat{a}}^2{e}^{i\delta t})$$

$$\frac{d}{dt}\left(\begin{array}{c}\hat{a}(t)\\ \hat{a}(t{)}^{†}\end{array}\right)$$

$$\frac{d^2}{d{t}^2}\hat{a}(t)+i\delta \frac{d}{dt}\hat{a}(t)-4{\kappa }^2\hat{a}(t)=0.$$

$$\hat{a}(\tau )=u_1(\tau )\hat{a}+u_2(\tau ){\hat{a}}^{†},$$

$$u_1(\tau )$$

$$u_2(\tau )$$

$$P[{\alpha }^{\prime };\tau ]=\frac{1}{{\pi }^2}\int d^2\beta e^{{\beta }^{\prime \ast }{\alpha }^{\prime }-{\beta }^{\prime }{\alpha }^{\prime \ast }}⟨:\hat{D}({\beta }^{\prime };\tau ):⟩.$$

$$⟨:\hat{D}(\beta ;\tau ):⟩\equiv \Phi (\beta ;\tau )=\exp \left[{\left(\begin{array}{c}{\beta }_r\\ {\beta }_i\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}\frac{1}{2}(1-|u_2^{\ast }-u_1{|}^2)& \mathrm{Im}[u_1{u}_2]\\ \mathrm{Im}[u_1{u}_2]& \frac{1}{2}(1-|u_2^{\ast }+u_1{|}^2)\end{array}\right)\left(\begin{array}{c}{\beta }_r\\ {\beta }_i\end{array}\right)\right].$$

$$\mathbf{S}=({\mathbf{\chi }}_{+},{\mathbf{\chi }}_{-})=\left(\begin{array}{cc}\frac{-b_{+}}{\sqrt{1+|b_{+}{|}^2}}& \frac{-b_{-}}{\sqrt{1+|b_{-}{|}^2}}\\ \frac{1}{\sqrt{1+|b_{+}{|}^2}}& \frac{1}{\sqrt{1+|b_{-}{|}^2}}\end{array}\right),$$

$$\Phi (\gamma ;\tau )=\exp \left[{\left(\begin{array}{c}{\gamma }_r\\ {\gamma }_i\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{cc}-c_{+}& 0\\ 0& -c_{-}\end{array}\right)\left(\begin{array}{c}{\gamma }_r\\ {\gamma }_i\end{array}\right)\right],$$

$$\text{with }{c}_{\pm }=-\frac{1}{2}(1-|u_2^{\ast }-u_1{|}^2)\pm 2\frac{|u_1{u}_2|}{1+b_{\pm }^2}.$$

$$P[\alpha ;\tau ]=$$

$$A_r$$

$$A_i$$

$$P_{\Omega }[\alpha ;\tau ]=$$

$$T(y,g)$$

$$P_{\Omega }[\alpha ;\tau ]=w^{2}T(w{A}_r,w^2{c}_{+})T(-w{A}_i,w^2{c}_{-}).$$

$$[\hat{a}(\tau ),\hat{a}(\tau +\Delta \tau )]$$

Full text

Quantum Optics

I Introduction

Nonclassical effects, such as photon antibunching KDM77 , squeezing W83 ; Slusher ; Wu ; Va08 ; Va16 , and entanglement EPR35 ; S35 , have been known for many decades. Their verification, classification, quantification, and application remain challenging tasks of modern quantum optics. For the distinction of genuine quantum interferences from classical optical effects, two major techniques have been established. The first one is the application of various types of nonclassicality criteria based on observable quantities. The second one is the investigation of different kinds of phase-space distributions. The implementation of each of these techniques brings along its own characteristic advantages and challenges.

A variety of nonclassicality criteria are suitable for different applications, depending on the quantum system or effect under consideration. Some of them consist of an infinite hierarchy of necessary and sufficient nonclassicality probes. Examples are criteria based on characteristic functions, moments, and their combination RiVo02 ; ShVo05 ; RSAMKHV15 . Yet, the full characterization of quantum effects requires the study of all orders of these hierarchies. Even though this is impossible in general, these methods provide a plethora of sufficient nonclassicality conditions to successfully identify various types of quantum effects; cf. Ref. Miranowic10 for an overview. However, low-order criteria may fail to uncover the nonclassical character of particular quantum states.

In such cases, the investigation of phase-space distributions may be advantageous. Prominent examples are the Husimi $Q$ function H40 , the Wigner function W32 , the Glauber-Sudarshan $P$ function S63 ; G63 , their unification in terms of $s$-parameterized quasiprobability distributions CG69 , and the general distributions introduced by Agarwal and Wolf AgWo . Nonclassicality is commonly defined via comparison of such quasiprobabilities with their classical counterparts. The widely accepted definition of nonclassicality by Titulaer and Glauber relies on the $P$ function TG65 ; M86 : Whenever $P$ cannot be interpreted as a classical probability density, i.e., when it contains negativities, the state is referred to as nonclassical. Hence, this very quasiprobability will be our benchmark to identify quantum effects. However, the study of the $P$ function can become a cumbersome task, as it is highly singular for many quantum states S16 . In general, an experimental reconstruction of the $P$ function is only possible if a proper regularization procedure is introduced KV10 . Based on such a technique, one can also implement a direct sampling approach for the regularized $P$ functions ASVKMH15 , which yields direct access to the full information on general quantum states.

The $P$ function itself represents the full information on the quantum state $\hat{\rho}$ of a radiation mode at an arbitrary but fixed time in a diagonal representation,

[TABLE]

with coherent states $|\alpha\rangle$ with complex amplitudes $\alpha$. Yet, correlations between multiple points in time (i.e., temporal correlations) play a fundamental role in quantum optics and quantum information theory. For instance, the early demonstration of the quantum nature of light via photon antibunching was based on the detection of two-time intensity correlation properties KDM77 . Also, the photocounting theory depends on multitime correlation functions of the light field to be measured KK64 . It is noteworthy that the corresponding field correlations are normal and time ordered. In particular, this is relevant when the interaction dynamics is described by an explicitly time-dependent Hamiltonian or when the radiation field is emitted by atomic sources. Recently, the Keldysh-ordered full counting statistics has been studied in the context of negativities of quasiprobabilities HC16 ; H17 .

Early studies of parametric processes in terms of time-dependent correlation functions were reported by Mollow Mollow73 . Temporal correlations have also been studied in the quantum dynamical theory of the fluctuations in a degenerate optical parametric oscillator PW94 . More recently, the multimode parametric dynamics has been further investigated by including time-ordering effects CBMS13 . For an explicitly time-dependent parametric interaction, temporal quantum correlations have also been considered in terms of characteristic functions KSV16 .

Another field where time ordering is important is the passive filtering of light emitted from atomic sources. A general theory for the effects of passive optical systems on quantum light was developed in Ref. KVW87 . In this context, a careful treatment of non-equal-time commutators is crucial. Although such commutation rules are not explicitly known in general, their effects can be handled in a closed form; for details see Chap. 2.7 in Ref. VW06 . Yet another example stems from the spectral filtering of quantum light from atomic sources for which spectral squeezing KVW86 and spectral intensity correlations Knoell86 ; Cresser87 of the atomic resonance fluorescence have been studied to some extent.

Temporal correlations are also considered in the framework of the Leggett-Garg inequalities LG85 ; ELN14 —sometimes referred to as temporal Bell inequalities. Recently, the latter were extended to continuous-variable systems placed in a squeezed state MV16 . It turns out that the application areas of temporal correlation properties of radiation fields are a wide-ranging field of research KK64 ; Mollow73 ; PW94 ; CBMS13 ; KSV16 ; KVW87 ; VW06 ; KVW86 ; Knoell86 ; Cresser87 ; LG85 ; ELN14 ; MV16 ; BKetal15 ; arxNMR ; cpetal15 . Hence, a complete characterization of such correlations is a subject of broad interest.

Also, for applications of time-dependent quantum correlations of light in quantum technology, a full characterization of such complex quantum effects is indispensable. For this purpose, a space-time-dependent phase-space representation has been introduced by generalizing the Glauber-Sudarshan $P$ function to a space-time-dependent $P$ functional V08 . This functional renders it possible to formulate an infinite set of nonclassicality conditions in terms of normal- and time-ordered field correlation functions, which are accessible by homodyne correlation measurements SV06 . However, such a verification of quantum correlations is hardly used, as it requires the detection of a manifold of correlation functions.

In conclusion, a direct study of the $P$ functional would be favorable and would lead to a deep and general understanding of temporal quantum correlations. However, even for a single time, the $P$ functional can be highly singular. Even more severely, the singularities of the multitime $P$ functional are even not clearly understood yet. Although it is well known that the singularities of the Glauber-Sudarshan $P$ function are caused by the normal-ordering prescription, it is—to our best knowledge—presently unknown whether or not the time-ordering prescription, occurring in the $P$ functional, can give rise to singularities stronger than those of a multimode (but single-time) $P$ function. In our recent contribution KSV16 , we formulated nonclassicality tests to uncover time-dependent quantum effects. This method was based on the characteristic function, i.e., the Fourier transform of the $P$ functional, which is indeed a regular function. Still, until now, a proper regularization procedure for the $P$ functional itself has not been established for the multitime case. Such a method, however, would be important for a comprehensive understanding and potential applications of general quantum correlations of radiation fields. Here, it is also noteworthy that the general quantum correlations under study even include entanglement as a subset SV09 .

In the present paper, we develop a rigorous formalism that describes nonclassical multitime correlations in terms of smooth nonclassicality quasiprobabilities. The regularity of our phase-space function applies to any evolution of the optical system. This enables us to verify quantum effects through the negativity of our quasiprobability representation for quantum states of general radiation fields. For example, we apply our method to an explicitly time-dependent parametric process and to a nonlinear dynamics of the quantized motion of a trapped ion. It is shown that the multitime scenario implies singularities of the $P$ functional stronger than those of the single-time multimode $P$ function. Even those singularities are suppressed by our regularization technique. Eventually, we formulate the measurement theory for the direct sampling of the regularized $P$ functional, which allows for an experimental visualization of multitime quantum correlations of light in phase space.

The paper is structured as follows. Section II recapitulates the concept of single- and multitime nonclassicality. In Sec. III, we formulate the regularization procedure of the multitime-dependent $P$ functional. Afterwards, in Sec. IV we apply the introduced techniques to an optical parametric process with frequency mismatch. The nonlinear evolution of a laser-driven trapped ion is analyzed in Sec V, which includes more complex time-dependent commutation rules and strongly enhanced singularities of the $P$ functional. In Sec. VI, we propose an experimental scheme for efficiently measuring the quantities under study. Finally, a summary and some conclusions follow in Sec. VII.

II Nonclassicality

II.1 Single-time nonclassicality filters

In the single-time scenario, the $P$ function can be used to express any quantum state as a formal mixture of coherent states G63 ; S63 ; TG65 ; M86 ; VW06 ; cf. Eq. (1). However, it is due to the singular behavior of $P$, which occurs for many quantum states, such as Fock or squeezed states, that a direct experimental access to this distribution is impossible. The singularity is equivalent to an unbounded characteristic function $\Phi$, being the Fourier transform of $P$.

To resolve the issue of singularities and to identify the nonclassicality of quantum states in experiments through negativities of properly defined quasiprobabilities, a regularization procedure was introduced KV10 . The resulting regularized $P$ function is consequently defined via the Fourier transform $\mathcal{F}$ or convolution $\ast$,

[TABLE]

with the filter function $\Omega_{w}$ or its inverse Fourier transform $\tilde{\Omega}_{w}$ and both depending on a width parameter $w>0$. This, in general, non-Gaussian filter function needs to satisfy several conditions KV10 :

$\Omega_{w}(\beta)\Phi(\beta)$ is rapidly decaying for all (finite) values $w$. This is necessary to assure that the regularized $P$ function exists and is even smooth for all states and for all filter widths (cf. also Ref. ASV13 , Appendix A). 2. 2.

The Fourier transform $\tilde{\Omega}_{w}$ is a probability density. Especially, its nonnegativity is important, as we want to visualize the negativities of the original $P$ function, which define the nonclassicality. Thus, the filter must not contribute negativities. 3. 3.

The limit $\lim_{w\rightarrow\infty}\Omega_{w}(\beta)=1$ assures that the original $P$ function is recovered for $w\to\infty$.

The regularized $N$-mode and single-time $P$ function is obtained from the generalization ASV13

[TABLE]

where $\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathbb{C}^{N}$. One possibility for constructing a multimode filter is a product of single-mode filters, $\tilde{\Omega}_{w}(\boldsymbol{\alpha})=\prod_{j=1}^{N}\tilde{\Omega}_{w}(\alpha_{j})$. Note that the multimode characteristic function $\Phi$ is, in general, unbounded, $\sup_{\boldsymbol{\beta}\in\mathbb{C}^{N}}|\Phi(\boldsymbol{\beta})|=\infty$, but can be bounded through a diverging function,

[TABLE]

which can be easily derived: Using the definition of the multimode characteristic function and the displacement operators $\hat{D}(\beta_{j})=\exp[\beta_{j}\hat{a}_{j}^{\dagger}-\beta_{j}^{\ast}\hat{a}_{j}]$ with $\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{N})^{\rm T}$, one obtains

[TABLE]

where $\hat{a}_{l}$ labels the annihilation operator of the $l$th radiation mode and $:\dots:$ connotes the normal ordering prescription. That is, all creation operators are placed to the left of the annihilation operators without making use of the bosonic commutation relations. As $[\hat{a}_{i},\hat{a}_{j}]=0$ and $[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\delta_{ij}$ hold ($\delta$ is the Kronecker symbol), one can use the standard (i.e., first-order) Baker-Campbell-Hausdorff (BCH) formula to obtain

[TABLE]

As for the unitary displacement operators $\hat{D}(\beta_{j})$ it holds that $\|\prod_{j=1}^{N}\hat{D}(\beta_{j})\|\leq 1$, and one readily verifies Eq. (4). This estimation also holds true for any time evolution of a quantum state $\hat{\rho}(t)$ for a single time $t$. Hence, one finds that the slope of the characteristic function of the $N$-mode $P$ function is bounded by an inverse Gaussian factor, which also bounds the singularities of the multimode $P$ function S16 . However, the regularization procedure so far recapitulated is restricted to single-time properties of a quantum system.

II.2 Multitime $P$ functional

In the more general multitime scenario, the situation is very different. It is nontrivial to give a similar expansion as in Eq. (1), because it is a cumbersome task to define the corresponding multitime density matrix APT14 . Thus, one needs a generalized, multitime-dependent version of the $P$ function V08 . We will discuss this concept in the continuation of this section. The resulting $P$ functional is formulated by using normal- and time-ordered expressions which are accessible in quantum correlation measurements VW06 ; SV06 . They also occur in the photocounting theory KK64 whenever source fields play a significant role in the description of a quantum state of light.

Let us consider an observable $\hat{O}[\{\hat{a}(t_{i})\}_{i=1}^{k}]$, which depends on the bosonic annihilation operators $\hat{a}(t_{i})$ and creation operators $\hat{a}(t_{i})^{\dagger}$ (not explicitly written as an argument in $\hat{O}$) at arbitrarily chosen points in time, $t_{1}\leq\dots\leq t_{k}$. Throughout this work, $k$ denotes the number of different points in time. For simplicity, we treat the case of a single spatial-frequency optical mode but $k$ nonmonochromatic (temporal) modes. The extension to $N$ spatial-frequency modes is straightforward, via $\hat{a}^{({\dagger})}\rightarrow(\hat{a}^{({\dagger})}_{1},\dots,\hat{a}^{({\dagger})}_{N})^{\rm T}$, as exemplified in the previous subsection for a single time. The definition of the singular $P$ functional was introduced in terms of space-time-dependent field operators in Ref. V08 , which already includes the most general scenario. It is a function of the $k$ coherent amplitudes at the considered $k$ points in time, $P[\alpha_{1},\dots,\alpha_{k};t_{1},\dots,t_{k}]$.

Using this $P$ functional, a multitime-dependent expectation value of the given observable $\hat{O}$ may be written as

[TABLE]

where we omitted the dependence on the complex conjugated variables and operators. The symbol $\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\dots\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}=\mathcal{T}:\dots:$ represents the normal ($:\dots:$)- and time ($\mathcal{T}$)-ordering prescription. Namely, the operators in Eq. (7) have to be normal ordered—creation operators to the left of annihilation operators—and then time ordering is performed, i.e., the time-dependent creation (annihilation) operators are sorted with increasing (decreasing) time arguments from left to right VW06 . Note that from the theory of photoelectric detection of light it is well known that observable correlation functions are subjected to normal and time ordering. An example of such a function is the second-order intensity correlation function, $g^{(2)}$, which corresponds to the expectation value of the operator ${\hat{O}\sim\hat{a}^{\dagger}(t)\hat{a}^{\dagger}(t+\Delta t)\hat{a}(t+\Delta t)\hat{a}(t)}$. This quantity is essential for the verification of photon antibunching KDM77 ; P82 .

From the general structure of Eq. (7), we can observe that the $P$ functional has formally the meaning of a joint probability distribution of the coherent amplitudes $\alpha_{i}\equiv\alpha(t_{i})$ at $k$ points in time. We use the term “formally” here, as $P$, in general, does not fulfill all the properties of a probability density in the sense of classical stochastics. This means that the $P$ functional is a joint quasiprobability, defined as the quantum expectation value V08

[TABLE]

where $\hat{\delta}$ denotes the operator-valued $\delta$ distribution. To clarify the terms, let us stress the following: The notion $P$ function is used for characterizing the quantum state at a single (arbitrary but fixed) time. The notion $P$ functional, on the other hand, applies when field amplitudes including their time dependencies are relevant. This also means that in the case $k=1$, the time-ordering prescription becomes meaningless, $\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\cdots\begin{smallmatrix}\circ\\ \circ\end{smallmatrix}\rightarrow{:}\,\cdots\,{:}$. In this scenario, the Glauber-Sudarshan $P$ function in Eq. (1) is recovered.

The classicality (nonnegativity) of the multitime functional, (8), leads to a hierarchy of classical inequalities in term of moments. Their violation certifies general quantum correlations of light V08 . Special cases were also studied in Ref. O07 for characterizing two-photon quantum interferences.

III The filtered $P$ functional

Starting from the definition (8), one can always express the $P$ functional through its characteristic function $\Phi$,

[TABLE]

where $\mathcal{F}_{\{\beta_{i}\}_{i=1}^{k}}=\mathcal{F}_{\beta_{1}}\cdots\mathcal{F}_{\beta_{k}}$ is a product of Fourier transforms for the $k$ different degrees of freedom. Let us recall that the (time-dependent) operator-valued $\delta$ distribution is defined as the Fourier transform of the (time-dependent) displacement operator, ${\hat{\delta}(\hat{a}(t)-\alpha)=\mathcal{F}_{\beta}[\hat{D}(\beta;t)](\alpha)}$, with $\hat{D}(\beta;t)=\exp\left[\beta\hat{a}(t)^{\dagger}-\beta^{\ast}\hat{a}(t)\right]$. Hence, one readily gets that

[TABLE]

which is the multitime-dependent characteristic function (MTCF) of the $P$ functional for ${t_{1}\leq\dots\leq t_{k}}$ KSV16 . The operator product for arbitrary $\hat{A}(t)$ is defined as $\prod_{i=1}^{k}\hat{A}(t_{i})=\hat{A}(t_{1})\dots\hat{A}(t_{k})$. Since the operators are in general not commuting, the operator products contain the time ordering from the first line in Eq. (III).

At first sight, the above expression, (III), resembles that for the multimode characteristic function [Eq. (5)]. However, there are two significant differences when considering the multitime scenario:

(i) One needs to consider non-equal-time commutators of the field operators, which are, in general, nonvanishing or not even proportional to unity VW06 ; C83 ; C84 ; KVW87 . This means that $[\hat{a}(t),\hat{a}(t^{\prime})^{\dagger}]\not\propto\hat{1}$ for $t\neq t^{\prime}$, which is a crucial point, as such commutators do not necessarily commute with other operators. Hence, the standard BCH formula is insufficient and higher-order terms need to be taken into account. Their calculation and the related convergence considerations are complex problems BC04 , which complicates the issue of finding a bound of the absolute square of the MTCF. This result is a major difference compared to the multimode (single-time) case [Eq. (4)].

(ii) Resulting from the structure of the functional, (8), the time-ordering prescription (beside the normal-ordering prescription) has to be considered. As we saw in the derivation of Eq. (4), the factor $e^{|\beta|^{2}/2}$ for a single radiation mode is caused by the normal ordering. The question arises whether or not the time ordering itself does cause a stronger rising behavior of the MTCF.

Altogether, the MTCF may be a more strongly growing function of $\beta_{1},\dots,\beta_{k}$ compared with the single-time but multimode scenario. However, this asymptotic behavior is important as it could lead to stronger singularities of the corresponding $P$ functional compared with the multimode $P$ function. This problem, to our best knowledge, has not been studied yet. Our rigorous regularization procedure in the multitime scenario has to include this eventuality. This also means that our approach, to be formulated, needs to be applicable to any dynamics of a quantum optical system.

III.1 Universal multitime regularization

As emphasized above, the singularities of $P$ in Eq. (9) are caused by the fact that the characteristic function $\Phi(\{\beta_{i};t_{i}\}_{i=1}^{k})$ [Eq. (III)] is, in general, unbounded. Hence, the integrals of the Fourier transforms in Eq. (9) do not converge. When generalizing the approach in Eq. (2), our multitime filter $\Omega_{\boldsymbol{w}}$ needs to assure the fast decay of the filtered characteristic function,

[TABLE]

Consequently, the regularized $P$ functional in terms of the Fourier transform is defined as

[TABLE]

where $\boldsymbol{w}$ denotes a tuple of width parameters that is specified later.

In the following, let us formulate some simple observations. We consider a continuous function $\varphi(z)$ depending on a real-valued parameter $z$ that might diverge for $|z|\to\infty$. In addition, we employ the triangular function,

[TABLE]

It holds that the product $\varphi(z)\mathrm{tri}(z)$ is bounded and continuous since the product of two continuous functions is continuous, and since the triangular function has the compact support $[-1,1]$, it is identical to [math] for $|z|>1$. Further on, it is easy to check that the one-dimensional Fourier transform of $\mathrm{tri}(z)$ is a probability density. Rescaling the argument with $w>0$, $\mathrm{tri}(z/w)$, yields a rescaled probability density, with the support of $\mathrm{tri}(z/w)$ changing to the interval $[-w,w]$. In particular, for $w\to\infty$, this rescaled triangular function converges pointwise to the constant function $1$, for which the Fourier transform is a $\delta$ distribution.

Returning to our initial filtering problem and keeping those observations in mind, we define the filter function

[TABLE]

using different filter parameters for each time, $\boldsymbol{w}=(w_{1},\dots,w_{k})$, and ${\rm Re}[\beta]$ with ${\rm Im}[\beta]$ denoting the real and imaginary part of $\beta$, respectively. From our observation above it directly follows that requirements 1–3 (cf. Sec. II) for a filter are satisfied. In particular and due to its compact support, the filter, (14), suppresses any rising behavior of the MTCF. This also means that the filtered $P$ functional, (12), exists always as a smooth function comment .

For illustration, a plot of the filter, multiplied by a factor $\exp[{\rm Re}[\beta]^{m}]$ for $m=2,4,6$, is given in Fig. 1. In the case $m=2$, the factor resembles the growth factor of a single-time-characteristic function P86 . The parameters $m=4,6$ correspond to faster increments that might result from non-equal-time commutation relations of nonlinear interactions. Note that a filter of this type is needed for multitime-dependent phenomena, as the behavior of the MTCF cannot be estimated in general. The multitime commutation rules of the bosonic operators are unknown for general interaction problems; for details see Chap. 2.7 of VW06 . With our filters in Eq. (14), however, we ensure that the multitime nonclassicality quasiprobabilities are smooth functions and they can be directly sampled in experiments, which is shown in Sec. VI.

III.2 Discussion

Let us formulate some preliminary conclusions. The problem of regularizing the $P$ functional has been treated. Although the singularities of this multitime quasiprobability are unknown, we formulated a regularization approach via a compact filter function, (14), which applies to any nonlinear interaction dynamics of radiation fields. On this basis, the well-behaved quasiprobability, (12), exists and it is negative for a width $w$ if and only if the light field under study is nonclassical—including multitime quantum correlations. Our applied filter suppresses any rising behavior of the MTCF for $|\beta_{l}|\to\infty$. This also includes scenarios where the slope of the MTCF increases more rapidly than an inverse Gaussian, which might result from the dynamics of complex interactions.

Remarkably, a triangular filter has already been used for the single-mode and -time scenario KV10 . There, the compact support was considered a deficiency because parts of the characteristic function are multiplied by [math] and, thus, do not contribute to the filtered $P$ function. Here, however, this compactness serves as a beneficial resource which renders it possible to filter a multitime $P$ functional. Moreover, similarly to the multimode scenario ASV13 , we have a filter in a product form, (14). It is worth mentioning that we could also and equivalently employ any filter with a compact support, including radial symmetric filters and higher-order autocorrelation-function filters KV12 ; KV14 .

The basic definition of the multitime nonclassicality is the inability to interpret the singular $P$ functional, (8), as a joint probability distribution of a classical stochastic process V08 . Also in Ref. V08 , a hierarchy of quantum correlation conditions has been formulated on the basis of measurable space-time-dependent field correlation functions. An approach to formulating multitime nonclassicality tests on the basis of the MTCF was subsequently formulated as well KSV16 . Here, in contrast, we introduce regular nonclassicality quasiprobabilities that enable us to study quantum correlations between multiple points in time directly via the corresponding negativities of smooth phase-space distributions of quantum stochastic, optical processes.

IV Parametric Processes

Let us now apply our approach to the characterization of temporal quantum effects on a specific physical system. In our recent work KSV16 , we studied the parametric process based on the MTCF. This and related parametric interactions are a fundamental tool for generating nonclassical light in modern experimental quantum optics, e.g., squeezed light Slusher ; Wu ; Va08 ; Va16 . Note that various single-photon sources are based on parametric down-conversion HO87 ; CS12 ; CD06 ; PJ05 ; RS04 ; LH01 . Due to the resulting wide range of applications of this process, let us reconsider this system from the perspective of the technique derived here.

The effective Hamiltonian of the quantum system in the interaction picture is

[TABLE]

with a positive frequency mismatch $\delta=\omega_{p}-2\omega_{a}$ and with $\omega_{p}$ and $\omega_{a}$ being the pump and signal frequency, respectively. In this process, a strong (classically described) pump field creates pairs of (equal-frequency) signal photons. Due to the violation of multitime-dependent classical inequalities, we have already demonstrated the presence of two-time quantum correlations KSV16 for this process. We also stress that the dynamical behavior exhibits a nontrivial dependence on time, as we have, in general, a nonvanishing commutator $[\hat{H}_{\rm int}(t),\hat{H}_{\rm int}(t^{\prime})]\neq 0$ for different times $t$ and $t^{\prime}$. In this section, let us focus on the visualization of quantum correlations directly in terms of negativities of the regularized $P$ functional, which has not been considered before.

The coupled equations of motion of the signal field operators $\hat{a}$ and $\hat{a}^{\dagger}$ read as

[TABLE]

After decoupling, one obtains second-order equations of motion,

[TABLE]

The solution can be found via standard algebra,

[TABLE]

where we have defined the following dimensionless quantities: $\vartheta_{r}=\pi\sqrt{16-r^{2}}/4$ (representing the eigenfrequency), $r=\delta/\kappa$ (the coupling ratio), $\tau=2\kappa t/\pi$ (a time in “natural” units of the system), and the two functions

[TABLE]

The initial condition is $\hat{a}(\tau=0)\equiv\hat{a}$.

IV.1 Single-time scenario

To clarify the filtering procedure and to demonstrate the applicability of the filter to the single-time dynamics of the system, we first consider the single-time scenario. In this case, $k=1$, the $P$ functional in Eq. (9) is obtained from the inverse Fourier transform, (III):

[TABLE]

Since we treat the time evolution in the interaction picture, the impact of the free-field Hamiltonian $\hat{H}_{0}=\hbar\omega_{a}\hat{a}^{\dagger}\hat{a}$ is not directly included in the parameters $u_{l}$. However, it only acts as a classical rotation in phase space, $\hat{D}(\beta^{\prime};\tau)\mapsto\hat{D}(\beta^{\prime}e^{i\omega_{a}t};\tau)$, which can be ignored, as we can perform a transformation $\beta^{\prime}e^{i\omega_{a}t}=\beta$ (likewise, $e^{i\omega_{a}t}\alpha^{\prime}=\alpha$ in the original phase space). Using Eq. (III) for $k=1$, inserting the solution of the time evolution, (18), using the decomposition $\beta=\beta_{r}+i\beta_{i}$ ($\beta_{r},\beta_{i}\in\mathbb{R}$), and reordering the terms using the BCH formula, we get

[TABLE]

Here, we have supposed that the initial state is the vacuum state and omitted to write the explicit time dependence, $u_{l}=u_{l}(\tau)$.

To simplify the integration, we diagonalize the coefficient matrix in Eq. (21). The transformation matrix $\boldsymbol{S}$, containing the normalized eigenvectors $\boldsymbol{\chi}_{l}$ ($l=+,-$), reads

[TABLE]

with $b_{\pm}=(-{\rm Re}[u_{1}u_{2}]\pm|u_{1}u_{2}|)/{\rm Im}[u_{1}u_{2}]\in\mathbb{R}$. The obtained normal coordinates are described via the classical rotation $(\beta_{r},\beta_{i})^{\rm T}=\boldsymbol{S}(\gamma_{r},\gamma_{i})^{\rm T}$ and we find

[TABLE]

Including the free-field propagation, we further obtain

[TABLE]

with

[TABLE]

The application of our regularization procedure, (12), to a single point in time yields

[TABLE]

The integral can be simplified (cf. Appendix C in S16 ) via defining the function

[TABLE]

which relates to complex error functions, and we finally arrive at

[TABLE]

The temporal evolution in terms of $P_{\Omega}[\alpha;\tau]$ is shown in Fig. 2 for different times $\tau=2\kappa t/\pi$. The negativities clearly display the nonclassicality of the system in terms of a regular and time-dependent quasiprobability.

IV.2 Singularities due to explicit time dependence

Let us now extend our studies to the more general multitime case. As discussed earlier, the main features of multitime correlations are due to (i) non-equal-time commutation relations and (ii) the time-ordering prescription; see the beginning of Sec. III. The commutators for (i) can be straightforwardly computed for the system under study by using the exact solution, (18),

[TABLE]

with $u_{l}(\tau)$ given in Eq. (IV). The effect of the time-ordering prescription in Eqs. (8) and (III) can be analyzed in terms of the ratio

[TABLE]

which relates the time- and normal-ordered quantities to solely normal-ordered ones. Since all appearing commutators are multiples of the identity, we get

[TABLE]

where $\varphi_{\beta_{j}}=\arg\beta_{j}$ for $j=1,2$. Using the solutions, (IV), one finds for $|r|\leq 4$ and $\vartheta_{r}=\pi\sqrt{16-r^{2}}/4$ that

[TABLE]

In other words, the influence of the time-ordering prescription is an additional term proportional to $\exp\{\pm|\beta_{1}\beta_{2}|\}$, where the sign depends on the sine ($\sin$) term in Eq. (IV.2). As the hyperbolic sine ($\sinh$) is monotonically increasing, the strength of this factor increases with $\Delta\tau$. For the $P$ functional (i.e., performing a Fourier transformation), this factor increases (for $\exp\{+|\beta_{1}\beta_{2}|\}$) or decreases (for $\exp\{-|\beta_{1}\beta_{2}|\}$) the strength of the singularities.

Let us apply our filter procedure introduced in Sec. III.1. We use identical filter widths, $w_{1}=w_{2}=w$ [cf. Eqs. (12) and (14)]. The general procedure to compute $P_{\Omega}[\alpha_{1},\alpha_{2};\tau_{1},\tau_{1}]$ is a straightforward extension of the one presented in Sec. IV.1. After some algebra, we get the two-time regularized $P$ functional,

[TABLE]

with the definition of $T$ in Eq. (27). The coefficients $f_{mn}$ and $d_{mn}$ together with a proper rotation of phase space can be obtained numerically as described in details in Sec. IV.1 and generalized to four dimensions.

The two-time regularized $P$ functional, (33), is depicted in Fig. 3. The negativities directly reveal the quantum correlations of the system under study. Note that for any time pairings $(\tau_{1},\tau_{2})$ and cross sections other than those used in Fig. 3, negativities are revealed as well; cf. also Ref. KSV16 . Let us stress that the existence of negativities for a certain filter width is necessary and sufficient for the existence of quantum correlations within the singular $P$ functional.

V Trapped-ion dynamics

In the previous section, we discuss multitime effects for a time-dependent parametric oscillator. In this case, however, the rather simple commutator $[\hat{a}(t),\hat{a}(t^{\prime})]\propto\hat{1}$ for $t\neq t^{\prime}$ holds [Eq. (IV.2)], and, hence, the standard BCH formula is applicable. Furthermore, one is able to solve the equations of motion analytically. The question arises how such systems shall be examined if the exact dynamics is not explicitly given. In this section, we therefore study more general structures of the MTCF and the corresponding dynamics for which the commutators are not central—i.e., they do not commute with the field operators itself, $[[\hat{a}(t),\hat{a}(t^{\prime})],\hat{a}(t^{(\prime)})]\neq 0$.

Let us first rewrite the MTCF [Eq. (III)] for the two-time case $k=2$ and an input state $\hat{\rho}_{\rm in}$,

[TABLE]

Here we have used $\hat{\rho}_{\rm in}(t_{1},t_{0})=\hat{U}(t_{1},t_{0})\hat{\rho}_{\rm in}\hat{U}(t_{1},t_{0})^{\dagger}$ and the time-evolved displacement operator $\hat{D}(\beta_{2};t_{2},t_{1})=\hat{U}(t_{2},t_{1})^{\dagger}\hat{D}(\beta_{2})\hat{U}(t_{2},t_{1})$. This form of $\Phi$ reveals a major difficulty: $e^{\beta_{1}\hat{a}^{\dagger}}$ and $e^{-\beta_{1}^{\ast}\hat{a}}$ are unbounded operators PB15 , and they cannot be rewritten in a simple manner when the commutators are not central.

Let us therefore consider the limit $t_{1}\to t_{0}$ and $t_{0}\to 0$, i.e., $t_{1}=t_{0}=0$. As in general $[\hat{a}(0),\hat{a}(t_{2})]\neq 0$ (even $\not\propto\hat{1}$) holds true, this situation differs from the single-time case. Let us emphasize that the time-ordering prescription still applies. Setting $t_{2}\equiv t$, we arrive at

[TABLE]

Evaluations of this expression depend on the input state $\hat{\rho}_{\rm in}$ and on the dynamics under study. Here we focus on Fock states as input states, i.e., $\hat{\rho}_{\rm in}=|p\rangle\langle p|$. This yields

[TABLE]

which is obtained via expanding the exponential functions in power series and using the standard actions of $\hat{a}$ ($\hat{a}^{\dagger}$) on the Fock states $|p\rangle$ ($\langle p|$).

First, we study Eq. (36) for $p=0$, i.e., $\hat{\rho}_{\rm in}=|{\rm vac}\rangle\langle{\rm vac}|$. We obtain

[TABLE]

which is always bound by an inverse Gaussian factor. Furthermore, $\Phi_{\rm vac}(\beta_{1},\beta_{2};0,t)$ equals a single-time characteristic function, and the corresponding $P$ functional [Eq. (9)] will attain for any dynamics the form

[TABLE]

As the delta distribution is a nonnegative distribution, the two-time $P$ functional describes a nonclassical system if and only if the single-time $P$ function fails to be a classical probability distribution. In this scenario, there are no genuine temporal correlations. However, the situation is different for other input states, which can be observed for our second example, $p=1$. Because $e^{-\beta_{1}^{\ast}\hat{a}}|1\rangle=|1\rangle-\beta_{1}^{\ast}|0\rangle$, we get

[TABLE]

and, therefore, additional terms due to the inclusion of a second point in time. This holds even if we set the first time to be [math]. Note, a similar behavior can be observed for any $p\geq 1$.

For clarification, let us consider a realistic interaction Hamiltonian WV97 ,

[TABLE]

which results in a time evolution obeying a noncentral commutator algebra and whose time evolution is solved numerically. It describes a nonlinear vibrational dynamics of a laser-driven trapped ion. Therein, $\varepsilon$ is the effective two-photon coupling strength, and $\eta$ is the Lamb-Dicke parameter. The nonlinear operator function $\hat{f}_{3}(\hat{a}^{\dagger}\hat{a};\eta)$ of the vibrational number operator $\hat{a}^{\dagger}\hat{a}$ accounts for the recoil effects due to absorption and emission of laser photons by the trapped atom. It reads WV97

[TABLE]

The action on Fock states of the ion’s center-of-mass motion yields

[TABLE]

with $L_{n}^{(k)}(x)$ being the generalized Laguerre polynomials. Using the completeness relation of the Fock states, $\sum_{n=0}^{\infty}|n\rangle\langle n|=\hat{1}$, the Hamiltonian, (40), can be written in the Fock basis as

[TABLE]

with $g_{3}(n;\eta)=f_{3}(n;\eta)\sqrt{(n+1)(n+2)(n+3)}$.

The time evolution can be numerically solved via evaluating the time evolution operator

[TABLE]

in matrix representation, where a sufficiently high cutoff of the Fock space has to be chosen. Using the MTCF in Eq. (35), we can investigate the difference between the squared modulus of the MTCF and the inverse Gaussian bound,

[TABLE]

with the dimensionless system time $\tau\equiv\varepsilon t$. If this function $\Delta\Phi$ exceeds [math], the MTCF is—due to temporal correlations—differently bounded compared to the two-mode single-time characteristic function. This means that the temporal correlations increase the divergences of the $P$ functional to an extent which cannot occur for any two-mode correlations at equal time; cf. Eq. (4).

A visualization of (45) is given in Fig. 4 for a particular choice of the parameters. We used $p=3$ ($\hat{\rho}_{\rm in}=|3\rangle\langle 3|$) and numerically evaluated $\Delta\Phi$ in a 200-dimensional Fock space to ensure approximation errors of the order of those of the numerical arithmetic. As for several parameters $\Delta\Phi>0$ holds, the temporal correlations obviously exceed the slope one can maximally expect for an equal-time two-mode characteristic function. As discussed earlier, the strength of the excess depends on the chosen input state. This deviation from the inverse Gaussian bound of the MTCF could be even stronger for other dynamical systems. However, our filter approach, introduced in Sec. III.1, is suitable for regularization of the $P$ functional for any dynamics. Based on the strongly nonlinear trapped-ion interaction Hamiltonian and the discussion of the impact of the input state, we have demonstrated our approach’s requirement for regularization of multitime $P$ functionals.

VI Sampling of the $P$ functional

The first experimental reconstruction of a phase-insensitive and single-time filtered $P_{\Omega}$ function with negativities was performed for single-photon-added thermal states KVBZ11 —based on the measurement of quadratures in balanced homodyne detection. A direct sampling of a single-time $P_{\Omega}$ function for a squeezed state was then performed KVHS11 —including the formulation of suitable pattern functions for discrete phase measurements. More recently, a method was presented to sample the filtered $P$ function via continuous-in-phase measurement ASVKMH15 . In the following we study an optical measurement scheme for reconstructing the filtered $P$ functional that allows one to apply our technique in experiments. Additionally note that a corresponding measurement technique can also be provided for the motional quantum state of a trapped ion, which can be directly based on the motional-state reconstruction as proposed in Ref. Wa-VO95 .

The setup of our scheme is shown in Fig. 5; it correlates the radiation field at two times. According to the quantum theoretical model for photodetection VW06 ; KK64 ; G65 , the joint probability of four detectors is

[TABLE]

where $n_{i}$ denotes the number of photons at the $i$th detector, $\eta_{i}$ is the detector efficiency, and $\hat{n}_{i}(t_{i})$ is the photon number operator at time $t_{i}$ in the corresponding detector path. As we correlate two points in time in our setup, we set $t_{1}=t_{2}\equiv t$ and $t_{3}=t_{4}\equiv t+\Delta t$. The phases of the local oscillator modes can be controlled separately by the phase shifters $\Delta\varphi$ and $\Delta\varphi^{\prime}$. Note that the layout is scalable and could be further extended to more points in time.

In close analogy to the procedure in Ref. VG93 , the correlated difference statistics is found to be

[TABLE]

with the difference events $n_{1}-n_{2}=v$ and $n_{3}-n_{4}=v^{\prime}$, the phases $\varphi$ and $\varphi^{\prime}$, the common amplitude $R$ of the two local oscillators, and the quadrature operator

[TABLE]

To arrive at Eq. (VI), we replaced $n_{i}$ with continuous variables, which can be done in the strong local oscillator limit and for a sufficiently large number of events. The distribution, (VI), is obviously the quantum expectation value of the time- and normal-ordered product of two Gaussian distributions of the quadratures at $t$ and $t+\Delta t$. Via the two-dimensional Fourier transform of the measured difference statistics,

[TABLE]

one gets

[TABLE]

Here $\Phi$ is the desired two-time characteristic function of the $P$ functional. Hence, by adjusting the parameters and properly scaling the complex numbers $y$ and $y^{\prime}$, $\Phi$ can be directly sampled with our setup. If we use the representations $\beta_{1}=b_{1}e^{i\varphi}$ and $\beta_{2}=b_{2}e^{i\varphi^{\prime}}$ and identify $y\eta R\equiv b_{1}$ and $y^{\prime}\eta^{\prime}R\equiv b_{2}$, we can rewrite the previous results as

[TABLE]

The regularized $P$ functional can be reconstructed in the following way. First, we recall that [cf. Eq. (12) for $k=2$]

[TABLE]

Since we have used a product filter [cf. Eq. (14)], we can simply insert Eq. (51) and, following the procedure in Ref. KVHS11 , rewrite the previous formula as

[TABLE]

Here, the so-called patten function $f_{\Omega}$ is

[TABLE]

with $z\in\{v,v^{\prime}\}$, $\eta_{v}\equiv\eta$, and $\eta_{v^{\prime}}\equiv\eta^{\prime}$. Finally, the regularized $P$ functional can be sampled from $M$ measured quadrature data points $(v_{j},\varphi_{j},v_{j}^{\prime},\varphi_{j}^{\prime})_{j=1}^{M}$ in the two channels via its empirical estimate

[TABLE]

where the time dependences are included in the set of data, $v\equiv v(t)$ and $v^{\prime}\equiv v(t+\Delta t)$. For convenience, we used a radial symmetric filter, such that the filter functions $\Omega_{w}$ depend only on the $b_{i}$ and not on the phases $\varphi$ and $\varphi^{\prime}$ [cf. Eqs. (VI) and (VI)].

Hence, we have formulated the theory of our proposed measurement scheme in Fig. 5 that renders it possible to directly obtain the two-time regularized $P$ functional via the sampling of measured (correlated) quadrature data using pattern functions. We may stress again that this approach is scalable to an arbitrary number of points in time by employing multiple BHDs. In addition, some remarks concerning the sampling error estimation can be found in Ref. ASVKMH15 .

VII Summary and Conclusions

In summary, we have derived a method for visualizing general multitime quantum correlations in terms of regular phase-space quasiprobabilities. For this purpose, we significantly generalized the approach of filtering the Glauber-Sudarshan $P$ function for a single time to a regularization of the $P$ functional that correlates an arbitrary number of points in time. This generalizes the commonly accepted definition of nonclassicality by Titulaer and Glauber to capture multitime quantum phenomena in terms of regular and accessible phase-space distributions. Hence, our formulation of multitime and regular nonclassicality quasiprobabilities is not restricted to particular correlation functions or observables. Additionally, we have proved that our method is applicable to arbitrarily complex evolutions of light fields. In our general approach, nonclassical correlations, if present, are directly visualized by negativities of this regular quasiprobability density for properly chosen filter widths. Beyond previously studied methods, our treatment enables us to visualize general quantum correlations of radiation fields, including quantum entanglement as a subset.

We applied this technique to characterization of the temporal quantum features of a parametric oscillator with a frequency mismatch. We studied the impact of the time-dependent Hamiltonian on the dynamical properties of this system, including the multitime correlations. In particular, we uncovered quantum correlations via a negative and regular two-time quasiprobability description of this process. We also studied a strongly nonlinear dynamics of the motional quantum state of a laser-driven trapped ion. In this case, nontrivial time-dependent commutator rules become important, leading to singularities of the $P$ functional much stronger than those occurring for equal-time two-mode correlations. Our technique regularizes these unexpectedly strong singularities. Eventually, we presented an experimental setup—consisting of two correlated balanced homodyne detection layouts. Based on the derived pattern functions, this allows one to directly sample the regularized quasiprobabilities in experiments. Altogether, this yields a powerful tool for the characterization of general, time-dependent quantum correlations in phase space.

Acknowledgements.

We thank Regina Kruse for helpful comments. This work was supported by the Deutsche Forschungsgemeinschaft through SFB 652, project B12. J.S. and W.V. acknowledge funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 665148.