# Precessing Binary Black Holes as Better Dark Sirens

Qianyun Yun1,2,Wen-Biao Han${}^{3,1,4,5},$ Qian Hu6, Haiguang Xu2
1Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai,
200030, China
2School of Physics and Astronomy, Shanghai Jiao Tong University 800 Dongchuan
RD.,Minhang District, Shanghai, 200240, China
3Hangzhou Institute for Advanced Study, University of Chinese Academy of
Sciences, Hangzhou 310124, China
4 School of Astronomy and Space Science, University of Chinese Academy of
Sciences, Beijing, 100049, China
5 Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou),
University of Chinese Academy of Sciences, Beijing 100049, China
6 Institute for Gravitational Research, School of Physics and Astronomy,
University of Glasgow, Glasgow, G12 8QQ, UK

(Accepted 2023 August 14. Received 2023 August 12; in original form 2023 June
13)

###### Abstract

Gravitational waves (GWs) from binary black hole mergers provide unique
opportunities for cosmological inference such as standard sirens. However, the
accurate determination of the luminosity distance of the event is limited by
the correlation between the distance and the angle between the binary’s
orbital angular momentum and the observer’s line of sight. In the letter, we
investigate the effect of precession on the distance estimation of binary
black hole events for the third-generation (3G) GW detectors. We find that the
precession can enhance the precision of distance inference by one order of
magnitude compared to the scenario where precession is absent. The constraint
on the host galaxies can be improved due to the improved distance measurement,
therefore the Hubble constant can be measured with higher precision and
accuracy. These findings underscore the noteworthy impact of precession on the
precision of distance estimation for 3G ground-based GW detectors, which can
serve as highly accurate probes of the Universe.

###### keywords:

gravitational waves – cosmological parameters – black hole mergers.

††pubyear: 2023††pagerange: Precessing Binary Black Holes as Better Dark
Sirens–LABEL:lastpage

##  1 Introduction

In 2015, the Advanced Laser Interferometer Gravitational-wave Observatory
(LIGO; (Harry, 2010) made a groundbreaking discovery in astronomy by detecting
gravitational wave (GW) event GW150914 ( Aasi et al. 2015; Abbott et al.
2016a, b, 2019b); this was followed in 2017 by the second Observation Run (O2)
detections (Acernese et al., 2015) from a binary neutron star system (BNS),
GW170817 (see, e.g., Abbott et al. 2019b, a, 2017b), which ignited the
emergence of multimessenger astronomy (MMA). The first multimessenger
observations of this merger provided the first standard siren measurement of
the Hubble constant, $H_{0}=69_{-8}^{+17}\ {\rm kms^{-1}Mpc^{-1}}$ (see, e.g.,
Guidorzi et al. 2017; Fishbach et al. 2019; Hotokezaka et al. 2019).

The discovery of the GW170817 was a significant breakthrough in the field of
gravitational wave astronomy. When a gravitational wave event has an
electromagnetic (EM) counterpart, the source redshift can be determined by EM
observations and be used to infer cosmological parameters, which is referred
as a ‘bright siren’ (Abbott et al., 2021). However, for events without
confirmed EM counterparts, such as the merger of binary black holes (BBHs), it
is challenging to determine the host galaxies and their redshifts directly.
‘Dark siren’ has been proposed (Schutz 1986; Chen et al. 2018; Fishbach et al.
2019; Finke et al. 2021; Leandro et al. 2022; Muttoni et al. 2023; Gair et al.
2023) as a statistical framework to infer cosmological parameters without an
EM counterpart. In the case of dark sirens, redshift information can be
obtained through statistical analysis of the gravitational wave events in
association with galaxy catalogs (see, e.g., Holz & Hughes 2005; MacLeod &
Hogan 2008; Cutler & Holz 2009; Taylor & Gair 2012; Finke et al. 2021;
Mancarella et al. 2022; Branchesi et al. 2023; Abbott et al. 2023). However,
compared to bright sirens, the constraints on cosmological parameters obtained
from dark sirens are much less stringent due to the lack of redshift
information. The large uncertainty in GW localization makes it difficult to
identify the true host galaxy and determine the redshift (see, e.g., Fishbach
et al. 2019; Soares-Santos et al. 2019; Palmese et al. 2020).

Third-generation gravitational (3G) wave detectors (Punturo et al., 2010),
like the Einstein Telescope (see, e.g., Punturo et al. 2010; Sathyaprakash et
al. 2012; Maggiore et al. 2020) and Cosmic Explorer (see, e.g., Reitze et al.
2019; Hall et al. 2021; Evans et al. 2021), will be a major improvement over
the current ground-based detectors. They will have much better sensitivity and
extend their low-frequency limit from 10 - 20 Hz to 1 Hz, allowing for the
detection of $10^{5}-10^{6}$ GW signals from compact binary coalescences over
a few years of observations (Sathyaprakash et al., 2011). These detectors will
provide a wealth of data for astrophysics, cosmology, and fundamental physics,
and will allow us to explore the dark energy sector by studying the equation
of state of dark energy and modified gravity theories (see, e.g., Dalal et al.
2006; Cutler & Holz 2009; Zhao et al. 2011; Cai & Yang 2017; Cai et al. 2018;
Belgacem et al. 2018; Yang et al. 2019; Branchesi et al. 2023; Yang et al.
2023).

In order to fully exploit the potential of dark sirens in cosmology, it is
crucial to accurately infer the luminosity distance and redshift of the source
(Gair et al., 2023). However, distance-inclination angle degeneracy is a
barrier to precise estimation of the source luminosity distance ($d_{L}$). A
recent publication (Borhanian et al., 2020) has demonstrated that events with
a large signal-to-noise ratio (SNR) could effectively break the degeneracy
between $d_{L}$ and $\iota$. Moreover, other studies have shown that some
physical effects, such as eccentricity and higher modes can help break the
degeneracy in GW parameter estimation (see, e.g., Yang et al. 2022, 2023; Ng &
et al 2023). Raymond et al. (2009) has shown that precession helps break
degeneracies in source sky direction for a two LIGO-like detector network.
(Raymond et al., 2009; Green et al., 2021) also pointed out that precession
can improve uncertainties in distance estimates. In this letter, we present
evidence that the precession of the compact binaries plays a role in
determining the distance for dark sirens for the 3G detectors. We anticipate
that 3G detectors will be capable of limiting the number of host galaxies for
precessing binary systems to a smaller range, whereas for non-precessing
binaries, the number could potentially be 2-10 times greater.

The presence of precession sheds light on the formation of compact binaries.
BBHs can form via two main mechanisms: common envelope evolution (Kalogera,
2000; Gerosa et al., 2018; Steinle & Kesden, 2022) and dynamical capture
(Benacquista & Downing, 2013; Rodriguez et al., 2016). In the case of
dynamical capture, the black hole spins are randomly oriented (Benacquista &
Downing, 2013; Rodriguez et al., 2016). The presence of precession sheds light
on the formation of compact binaries and also offers an opportunity for
testing black holes (Zhang et al., 2022). Precessing spins are predicted for
BBHs formed in dense clusters, but can also occur in isolated sources due to
supernova kicks. Using precessing BBHs as dark sirens requires accounting for
theoretical uncertainties in their number and distribution. Their spin
precession, influenced by factors like tides, winds, and high-mass binary star
accretion, can be induced with >20 per cent accretion of their companion’s
envelope. Even 2 per cent accretion from a Wolf-Rayet star as the accretor can
cause significant precession (Steinle & Kesden, 2022). Additionally,
hierarchical mergers of stellar-mass black holes are expected to have spins
clustered around 0.7 (Gerosa & Fishbach, 2021).

The goal of this letter is to investigate if spin precession in a BBH system
can increase the precision of distance estimations for the parameter
estimation of BBH signals. This letter endeavors to scrutinize BBH sources
situated within a luminosity distance span of 3,000 to 40,000 Mpc (redshift
$z$ = 0.5 - 5). We use
$d_{\mathrm{L}}(z)=c(1+z)\int_{0}^{z}\frac{dz^{\prime}}{H\left(z^{\prime}\right)}$
to converse luminosity distance to redshift z. The $H_{0}$ is 67.7 km /
(Mpc$\cdot$ s). The $\Omega_{m}$ is 0.3085. The $\Omega_{\Lambda}$ is 0. The
$\Omega_{k}$ is -1. As illustrated in (Frieman et al., 2008), during epochs
marked by a redshift range exceeding 0.5, the Universe is pre-dominantly
governed by matter. In contrast, redshifts below approximately 0.5 signify the
dominance of dark energy. As a result, observing and analyzing BBH sources at
these particular redshifts provides us with a precious understanding of the
shift from a matter-dominated to a dark energy-dominated phase in the
Universe.

##  2 Methodology

We employ the IMRPhenomPv3 model to produce precessing waveforms. IMRPhenomPv3
model takes the effects of two independent spins in the precessing dynamics
into account (Khan et al., 2019). This waveform approximant is available in
LALSuite (LIGO Scientific Collaboration, 2020). To create the precessing
waveform, the first step of IMRPhenomPv3 is to express the two polarizations
of the gravitational wave ($h_{+}$ and $h_{\times}$) as a linear combination
of spherical harmonics with spin-weight -2 (see, e.g., Pan et al. 2014; Hannam
et al. 2014). The GW multipoles $(h_{l,m}$) are the coefficients of these
basis functions, calculated by Equation. 1. This is done in a frame aligned
with the total angular momentum, and the direction of propagation is
determined by the spherical polar coordinates $(\theta,\phi)$.

$h(t;\theta,\phi)=h_{+}-ih_{\times}=\sum_{\ell\geqslant
2,m}h_{\ell,m}(t)Y_{\ell,m}^{-2}(\theta,\varphi).$ (1)

The IMRPhenomPv3 expresses the waveforms of precessing BBHs $h^{\rm
prec}_{l,m}$ as a function of the multipoles from non-precessing BBHs $h^{\rm
non-prec}_{l,m}$. We only include the dominant $\ell=2$ modes, i.e., (Khan et
al., 2019)

$h_{2,m}^{\text{prec
}}(t)=e^{-im\alpha(t)}\sum_{\left|m^{\prime}\right|=2}e^{im^{\prime}\epsilon(t)}d_{m^{\prime},m}^{2}(-\beta(t))h_{2,m^{\prime}}^{\text{non-
prec }}(t).$ (2)

Fisher matrix formalism is the common tool for forecasting detector
capabilities, however, it could lose validity in precessing binaries due to
the unconstrained spin angles and ill-behaved matrices (Borhanian &
Sathyaprakash, 2022). Therefore, we utilize the Bayesian parameter estimation
package Bilby (Ashton et al., 2019) to estimate the source parameters and
compare the localization inference capabilities of the different spins. The
standard Gaussian noise likelihood function $L$ is used to analyze the strain
data $d_{k}$ based on the source parameters: $\theta$ (see, e.g., Ashton et
al. 2019; Romero-Shaw et al. 2020):

$\ln\mathcal{L}(d\mid\theta)=-\frac{1}{2}\sum_{k}\left\\{\frac{\left[d_{k}-\mu_{k}(\theta)\right]^{2}}{\sigma_{k}^{2}}+\ln\left(2\pi\sigma_{k}^{2}\right)\right\\}.$
(3)

In the present work, we employ IMRPhenomPv3 to generate the waveform
$\mu(\theta)$, where $\theta$ consists of 15 parameters,
$\theta=\\{\mathcal{M}_{c},\rm
q,d_{L},\theta_{JN},ra,dec,\psi,a_{1},a_{2},\theta_{1},\theta_{2},\phi_{JL},\phi_{12},t_{c},\phi_{c}$},
$\mathcal{M}_{c}$ and $q$ are the chirp mass and mass ratio of the two black
holes. $d_{L}$ is the luminosity distance. ‘ra’ and ‘dec’ describe the sky
position of the event and $\psi$ is the polarization angle. $\rm\theta_{JN}$
is the inclination angle at the reference frequency (which is set to 10 Hz in
this work) for a processing system. $a_{1},a_{2}$ are the dimensionless spin
magnitudes of two black holes. The four angles
$\rm\theta_{1},\theta_{2},\phi_{JL},\phi_{12}$ quantify the spin misalignment
of the binary, which drives the system to precess. $t_{c}$ is the merging time
and $\phi_{c}$ is the coalescence phase. We marginalize the parameters $t_{c}$
and $\phi_{c}$ in this letter.

Figure 1: The posterior distributions of the luminosity distance for different
injected distances. The red line in each subfigure represents injected
distance for each case. $\rm\theta_{JN}$ is $\pi/4$ for all cases. For the
precession case, $\rm\phi_{JL}$ and $\phi_{12}$ are both 0.7, respectively.
The box in the middle of the violin plots represents the upper and lower
quartiles of the posterior distributions.

We inject a signal into Gaussian noise generated with ET-D (Punturo et al.,
2010) and CE (Abbott et al., 2017a) design sensitivity. We assume two
detectors are located at Virgo and LIGO-Hanford sites, respectively. It is
important to note that this study specifically focuses on a detector network
composed of ET and CE. We employ a fast sampler Nessai (Williams et al., 2021)
to generate posterior samples. Due to the inefficiency of samplers in high-SNR
scenarios, where the SNR can reach values above 180 for CE detector, a
significant number of computational resources will be consumed to perform
complete parameter estimation. We fix the spin parameters $\rm
a_{1},a_{2},\theta_{1},\theta_{2},\phi_{JL}$ and $\phi_{12}$ to injected
values in parameter estimation. The corresponding conditional posterior could
be narrower than a fully marginalized posterior, but it suffices to
investigate the distance-inclination degeneracy. The rationale for this
approach is twofold: Firstly, the high SNR ratio allows us to tightly
constrain the spin parameters, whereas the third Observation Run’s signal-to-
noise ratio effectively constrained precession (Hannam et al., 2022).
Secondly, by assuming perfect knowledge of spins, we obtain conditional
posterior distributions instead of marginalized posterior distributions, which
may result in overconstrained error bars. Although specific numerical values
might be affected, the fundamental trend of breaking the distance-inclination
degeneracy due to spin effects remains valid.

We inject two types of sources with the same component masses into the
detectors, one with precession and one without, to evaluate the detectors’
ability to locate them and to test whether precession can alleviate the
degeneracy between luminosity distance and inclination angle. The prior
distribution for the luminosity distance is assumed to be uniform within the
range of 1000-21000 Mpc for low redshift cases (z = 0.5-2) and 21000-70000 Mpc
for high redshift cases (z = 2-5). We employ Hellinger distance to compute the
difference between two probability distributions of priors and the estimated
results. The Hellinger distance (Hellinger, 1909; Moore & Gerosa, 2021)
between two discrete probability distributions p and q is : $\rm
d_{H}^{2}=1-\sum_{i=1}^{N}\sqrt{p_{i}q_{i}}$. Hellinger distance is a measure
of the similarity between probability distributions. Precessing and non-
precessing systems are set at the same luminosity distance and we investigate
various distances. The parameters used for injection are as follows: the chirp
mass is 50; the mass ratio is 1; and $a_{1},a_{2}$ are randomly selected from
the range of 0 to 0.9.

Fig.1 displays the posterior distribution of luminosity distance for the
precessing and non-precessing sources at different distances. In Fig. 1, the
violin plots visually represent the distribution of luminosity for sources at
different luminosity distances. The central lines represent the box plot data,
displaying the 25th, 50th, and 75th percentiles. The thin lines within the
plot represent the 90 per cent confidence intervals. It is obvious that the 90
per cent confidence intervals of the distance estimation for processing BBHs
are narrower compared to that of non-precessing BBHs. As depicted in Fig. 1,
precession can greatly improve the precision of luminosity distance
estimation, especially at large distances. With precession, even for sources
at greater distances, the luminosity distance can be estimated precisely. In
contrast, the uncertainties of $\Delta d_{L}$ becomes extremely large in
systems without precession. Fig. 1 also show the Hellinger distances for
different situations.

In Fig. 2, we explore the differences in distance estimation between
precessing and non-precessing BBHs by injecting 200 sources. These injected
BBHs have randomly assigned parameters. Their luminosity distance ranges from
3000 to 48000 Mpc, $a_{1}$ and $a_{2}$ range from 0 to 0.9. The presence of
precession results in a reduction of the confidence interval for the estimated
distance ($\rm\Delta d_{L}$). The $\rm\Delta d_{L}/d_{L}$ of precessing BBHs
is consistently below 10 per cent, whereas the $\rm\Delta d_{L}/d_{L}$ of non-
precessing BBHs exceeds 40 per cent. The improvement in the uncertainty of
distance estimation for precessing BBHs is more prominent at shorter distances
and gradually diminishes as the redshift increases, in line with the expected
decrease in SNR ratio with greater distances.

Figure 2: The confidence interval of luminosity distance divided by the
luminosity distance ($\rm\Delta d_{L}/d_{L}$) for 200 random sources at
different luminosity distances for both precessing cases (red dots) and non-
precessing cases (blue dots). The other parameters are set randomly.

Figure 3:  The comparison of parameter estimation result of precessing BBHs
and non-precessing BBHs. The injected luminosity distance is 11800 Mpc and
37974 Mpc. $\rm\theta_{JN}$ is $\pi/4$ for both case. For the two preccesing
BBHs, their $\rm\phi_{JL}$ and $\phi_{12}$ are both 1.2, respectively.
Distinct and randomized sky locations are utilized for the upper and lower
plots. Figure 4:  The comparison between the estimation of $H_{0}$ for dark
sirens originating from precessing BBHs and those from non-precessing BBHs.

Fig.3 displays an example of the parameter estimation outcome for sources
placed at a luminosity distance of 11800 Mpc. The non-precessing system shows
strong degeneracy between inclination angle and distance, which leads to an
inaccurate estimation of distance. The precession significantly alleviates the
degeneracy between $\theta_{JN}$ and the luminosity distance. The estimated
distance of a precessing BBH is more precise compared to a non-precessing BBH.

Our findings highlight the role of precession in the accurate localization of
distant gravitational wave sources. These results emphasize the importance of
precessing source when using GW as dark sirens to measure the Hubble constant
($H_{0}$) in future 3G detectors.

##  3 Discussion and Conclusion

In order to utilize GW events as standard sirens for measuring the Hubble
constant, it is essential to measure distances with high precision for the
identification of the host galaxy. For the dark sirens, the strong degeneracy
between the luminosity distance and the inclination angle renders it difficult
to infer the distance precisely. We emphasize the crucial role of precession
in localizing gravitational wave events. By performing Bayesian parameter
estimation, we demonstrate that precession significantly narrows the
uncertainty of estimating luminosity distance by up to 10 times compared to
non-precessing cases, especially for distant sources. This is because the
precession could largely alleviate the degeneracy between distance and
inclination angle. With the posterior localization of GW sources, one can
obtain the upper and lower limits of the redshift (see, e.g., Chen et al.
2018; Fishbach et al. 2019), then one can filter out the candidate host
galaxies. Finally, the value of $H_{0}$ for cosmic expansion can be obtained
through the Bayesian inference (see, e.g., Zhu et al. 2022; Song et al. 2022).
The precessing BBHs can reduce the uncertainty of the luminosity distance
greatly compared to the non-precessing cases, then the number of the host
galaxy candidates will be reduced too as it is roughly proportional to the
uncertainty of distance (Fan et al., 2014). Therefore, the measured Hubble
constant will have smaller uncertainty. We make the assumption of galaxies
being uniformly distributed within the co-moving volume (Barausse, 2012; Wang
et al., 2022). This assumption enables us to estimate the number of galaxies
within the given volume. By utilizing the inferred sky positions and distance
inferences, we calculate the comoving volumes in which the potential host
galaxies may reside. The number density of galaxies in the Universe is 0.02
$\rm Mpc^{-3}$ (Wang et al., 2022). Multiplying this density number by the
volume allows us to derive the count of potential host galaxies. This
calculation provides an estimation of the number of galaxies that have the
potential to serve as hosts. The inclusion of precession in BBH systems can
lead to a substantial decrease in the potential number of host galaxies, in
comparison to non-precessing BBH systems, typically by a factor of 2-10. This
reduction becomes particularly prominent when the sources are located at
distances greater than 10,000 Mpc.Additionally, we performed a simple
evaluation to investigate how the precession-distance effect could influence
the estimation of $H_{0}$ for the binary black hole sources depicted in Fig.
2. Our analysis uses the redshift of the host galaxy by converting the
injected distance into a corresponding redshift value. Subsequently, we
incorporated the 90 per cent interval of the estimated luminosity distance and
applied the standard cosmological model to calculate the range of $H_{0}$. The
expression for $H(z)$ was given by
$H(z)=H_{0}\sqrt{\Omega_{\mathrm{m}}(1+z)^{3}+1-\Omega_{\mathrm{m}}}$, where
the equation of state (EoS) parameter of dark energy, $w$, was assumed to be
-1. The values of $\Omega_{m}$, $\Omega_{k}$, and $\Omega_{\Lambda}$ remain
consistent with those mentioned earlier.

The results of our calculations are illustrated in Fig. 4. The blue curve
represents $H_{0}$ inferred from the precessing sources as dark sirens, while
the green curve represents $H_{0}$ inferred from the non-precessing sources.
It is evident from the figure that the precessing BBH sources can provide a
more tightly constrained estimation of $H_{0}$ as dark sirens. The
distributions in Fig. 4 are normalized to have equal area. The areas are equal
to 1.

The results in this letter suggest that precessing GW sources as dark sirens
have immense potential for understanding the cosmic expansion history. By
studying the expansion history, we can gain insights into the behavior of dark
energy and modified gravity theories (Belgacem et al., 2018). Although our
research primarily focuses on the third generation of ground-based
gravitational wave detectors, the results can be extended to space-based
detectors as well (Lang & Hughes, 2006; Lang et al., 2011; Klein et al.,
2009). The use of precessing dark sirens could pave the way for more accurate
and precise measurements of the Hubble constant and a deeper understanding of
the universe’s evolution.

## Acknowledgements

This work was supported by the National Key R&D Program of China (Grant Nos.
2021YFC2203002), and the National Natural Science Foundation of China (Grant
No. 12173071). Wen-Biao Han was supported by the CAS Project for Young
Scientists in Basic Research (Grant No. YSBR-006).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the
corresponding author Wen-Biao Han.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Class. Quant. Grav., 32, 074001
  * Abbott et al. (2016a) Abbott B. P., et al., 2016a, Phys. Rev. X, 6, 041015
  * Abbott et al. (2016b) Abbott B. P., et al., 2016b, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, Class. Quant. Grav., 34, 044001
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2019a) Abbott B. P., et al., 2019a, Phys. Rev. X, 9, 011001
  * Abbott et al. (2019b) Abbott B. P., et al., 2019b, Astrophys. J., 875, 160
  * Abbott et al. (2021) Abbott B. P., et al., 2021, Astrophys. J., 909, 218
  * Abbott et al. (2023) Abbott R., et al., 2023, Astrophys. J., 949, 76
  * Acernese et al. (2015) Acernese F., et al., 2015, Class. Quant. Grav., 32, 024001
  * Ashton et al. (2019) Ashton G., et al., 2019, Astrophys. J. Suppl., 241, 27
  * Barausse (2012) Barausse E., 2012, Mon. Not. Roy. Astron. Soc., 423, 2533
  * Belgacem et al. (2018) Belgacem E., Dirian Y., Foffa S., Maggiore M., 2018, Phys. Rev. D, 97, 104066
  * Benacquista & Downing (2013) Benacquista M. J., Downing J. M. B., 2013, Living Rev. Rel., 16, 4
  * Borhanian & Sathyaprakash (2022) Borhanian S., Sathyaprakash B. S., 2022, Listening to the Universe with Next Generation Ground-Based Gravitational-Wave Detectors (arXiv:2202.11048)
  * Borhanian et al. (2020) Borhanian S., Dhani A., Gupta A., Arun K. G., Sathyaprakash B. S., 2020, Astrophys. J. Lett., 905, L28
  * Branchesi et al. (2023) Branchesi M., et al., 2023, JCAP, 07, 068
  * Cai & Yang (2017) Cai R.-G., Yang T., 2017, Phys. Rev. D, 95, 044024
  * Cai et al. (2018) Cai R.-G., Liu T.-B., Liu X.-W., Wang S.-J., Yang T., 2018, Phys. Rev. D, 97, 103005
  * Chen et al. (2018) Chen H.-Y., Fishbach M., Holz D. E., 2018, Nature, 562, 545
  * Cutler & Holz (2009) Cutler C., Holz D. E., 2009, Phys. Rev. D, 80, 104009
  * Dalal et al. (2006) Dalal N., Holz D. E., Hughes S. A., Jain B., 2006, Phys. Rev. D, 74, 063006
  * Evans et al. (2021) Evans M., et al., 2021, A Horizon Study for Cosmic Explorer: Science, Observatories, and Community (arXiv:2109.09882)
  * Fan et al. (2014) Fan X., Messenger C., Heng I. S., 2014, Astrophys. J., 795, 43
  * Finke et al. (2021) Finke A., Foffa S., Iacovelli F., Maggiore M., Mancarella M., 2021, JCAP, 08, 026
  * Fishbach et al. (2019) Fishbach M., et al., 2019, Astrophys. J. Lett., 871, L13
  * Frieman et al. (2008) Frieman J., Turner M., Huterer D., 2008, Ann. Rev. Astron. Astrophys., 46, 385
  * Gair et al. (2023) Gair J. R., et al., 2023, Astron. J., 166, 22
  * Gerosa & Fishbach (2021) Gerosa D., Fishbach M., 2021, Nature Astron., 5, 749
  * Gerosa et al. (2018) Gerosa D., Berti E., O’Shaughnessy R., Belczynski K., Kesden M., Wysocki D., Gladysz W., 2018, Phys. Rev. D, 98, 084036
  * Green et al. (2021) Green R., Hoy C., Fairhurst S., Hannam M., Pannarale F., Thomas C., 2021, Phys. Rev. D, 103, 124023
  * Guidorzi et al. (2017) Guidorzi C., et al., 2017, Astrophys. J. Lett., 851, L36
  * Hall et al. (2021) Hall E. D., et al., 2021, Phys. Rev. D, 103, 122004
  * Hannam et al. (2014) Hannam M., Schmidt P., Bohé A., Haegel L., Husa S., Ohme F., Pratten G., Pürrer M., 2014, Phys. Rev. Lett., 113, 151101
  * Hannam et al. (2022) Hannam M., et al., 2022, Nature, 610, 652
  * Harry (2010) Harry G. M., 2010, Class. Quant. Grav., 27, 084006
  * Hellinger (1909) Hellinger E., 1909, Journal für die reine und angewandte Mathematik, 1909, 210
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, Astrophys. J., 629, 15
  * Hotokezaka et al. (2019) Hotokezaka K., Nakar E., Gottlieb O., Nissanke S., Masuda K., Hallinan G., Mooley K. P., Deller A. T., 2019, Nature Astron., 3, 940
  * Kalogera (2000) Kalogera V., 2000, Astrophys. J., 541, 319
  * Khan et al. (2019) Khan S., Chatziioannou K., Hannam M., Ohme F., 2019, Phys. Rev. D, 100, 024059
  * Klein et al. (2009) Klein A., Jetzer P., Sereno M., 2009, Phys. Rev. D, 80, 064027
  * LIGO Scientific Collaboration (2020) LIGO Scientific Collaboration 2020, LALSuite: LIGO Scientific Collaboration Algorithm Library Suite, Astrophysics Source Code Library, record ascl:2012.021 (ascl:2012.021)
  * Lang & Hughes (2006) Lang R. N., Hughes S. A., 2006, Phys. Rev. D, 74, 122001
  * Lang et al. (2011) Lang R. N., Hughes S. A., Cornish N. J., 2011, Phys. Rev. D, 84, 022002
  * Leandro et al. (2022) Leandro H., Marra V., Sturani R., 2022, Phys. Rev. D, 105, 023523
  * MacLeod & Hogan (2008) MacLeod C. L., Hogan C. J., 2008, Phys. Rev. D, 77, 043512
  * Maggiore et al. (2020) Maggiore M., et al., 2020, JCAP, 03, 050
  * Mancarella et al. (2022) Mancarella M., Borghi N., Foffa S., Genoud-Prachex E., Iacovelli F., Maggiore M., Moresco M., Schulz M., 2022, PoS, ICHEP2022, 127
  * Moore & Gerosa (2021) Moore C. J., Gerosa D., 2021, Phys. Rev. D, 104, 083008
  * Muttoni et al. (2023) Muttoni N., Laghi D., Tamanini N., Marsat S., Izquierdo-Villalba D., 2023, Dark siren cosmology with binary black holes in the era of third-generation gravitational wave detectors (arXiv:2303.10693)
  * Ng & et al (2023) Ng K. K. Y., et al 2023, Physical Review D, 107, 024041
  * Palmese et al. (2020) Palmese A., et al., 2020, Astrophys. J. Lett., 900, L33
  * Pan et al. (2014) Pan Y., Buonanno A., Taracchini A., Kidder L. E., Mroué A. H., Pfeiffer H. P., Scheel M. A., Szilágyi B., 2014, Phys. Rev. D, 89, 084006
  * Punturo et al. (2010) Punturo M., et al., 2010, Class. Quant. Grav., 27, 194002
  * Raymond et al. (2009) Raymond V., van der Sluys M. V., Mandel I., Kalogera V., Rover C., Christensen N., 2009, Class. Quant. Grav., 26, 114007
  * Reitze et al. (2019) Reitze D., et al., 2019, Bull. Am. Astron. Soc., 51, 035
  * Rodriguez et al. (2016) Rodriguez C. L., Zevin M., Pankow C., Kalogera V., Rasio F. A., 2016, Astrophys. J. Lett., 832, L2
  * Romero-Shaw et al. (2020) Romero-Shaw I. M., et al., 2020, Mon. Not. Roy. Astron. Soc., 499, 3295
  * Sathyaprakash et al. (2011) Sathyaprakash B., et al., 2011, in 46th Rencontres de Moriond on Gravitational Waves and Experimental Gravity. pp 127–136 (arXiv:1108.1423)
  * Sathyaprakash et al. (2012) Sathyaprakash B., et al., 2012, Class. Quant. Grav., 29, 124013
  * Schutz (1986) Schutz B. F., 1986, Nature, 323, 310
  * Soares-Santos et al. (2019) Soares-Santos M., et al., 2019, Astrophys. J. Lett., 876, L7
  * Song et al. (2022) Song J.-Y., Wang L.-F., Li Y., Zhao Z.-W., Zhang J.-F., Zhao W., Zhang X., 2022, Synergy between CSST galaxy survey and gravitational-wave observation: Inferring the Hubble constant from dark standard sirens (arXiv:2212.00531)
  * Steinle & Kesden (2022) Steinle N., Kesden M., 2022, Phys. Rev. D, 106, 063028
  * Taylor & Gair (2012) Taylor S. R., Gair J. R., 2012, Phys. Rev. D, 86, 023502
  * Wang et al. (2022) Wang R., Ruan W.-H., Yang Q., Guo Z.-K., Cai R.-G., Hu B., 2022, Natl. Sci. Rev., 9, nwab054
  * Williams et al. (2021) Williams M. J., Veitch J., Messenger C., 2021, Phys. Rev. D, 103, 103006
  * Yang et al. (2019) Yang T., Holanda R. F. L., Hu B., 2019, Astropart. Phys., 108, 57
  * Yang et al. (2022) Yang T., Cai R.-G., Cao Z., Lee H. M., 2022, Phys. Rev. Lett., 129, 191102
  * Yang et al. (2023) Yang T., Cai R.-G., Cao Z., Lee H. M., 2023, Phys. Rev. D, 107, 043539
  * Zhang et al. (2022) Zhang C., Han W.-B., Yang S.-C., 2022, Mon. Not. Roy. Astron. Soc., 516, L107
  * Zhao et al. (2011) Zhao W., Van Den Broeck C., Baskaran D., Li T. G. F., 2011, Phys. Rev. D, 83, 023005
  * Zhu et al. (2022) Zhu L.-G., Hu Y.-M., Wang H.-T., Zhang J.-d., Li X.-D., Hendry M., Mei J., 2022, Phys. Rev. Res., 4, 013247


Reproducing results in publications by distributing publicly available source
code is becoming ever more popular.  Given the difficulty of reproducing   experiments, there have been significant efforts in reducing the variance of
these results.  As in any science, the ability to consistently reproduce
results effectively strengthens the underlying hypothesis of the work, and
thus, should be regarded as important as the novel aspect of the research
itself.  The contribution of this work is a framework that is able to
reproduce consistent results and provides a means of easily creating,
training, and evaluating   models.

Consistently reproducing results is a fundamental criteria of the scientific
method, without which, a hypothesis may be weakened or even invalidated. This
is becoming even more necessary, because a growing number of publications are
inflated by false positives, to the point that they are labeled with the
pejorative term “p-hacking”, the intentional or unintentional act to bias
results in favor of publication for peer
acceptance <cit.>; efforts in this direction
include introducing new statistical methods to detect false findings

The inability to reproduce results has been referred to as the “replication
crisis” <cit.>.  The problem of
reproducibility in results is becoming more acknowledged as a serious issue in
the  community with efforts to understand and overcome the
challenge <cit.>. Not only has
the community addressed the issue in the literature, it has endeavored to
assess if experiments are reproducible and provide recommendations to remedy
the problem where reproducibility is lacking.  An example of this effort
includes reporting on experimental methodology, implementation, analysis and
conclusions in the .

To address these issues, we present , a  framework for  research
by and for the academic research community.  Not only does the framework
address issues of reproducibility, it also is designed to easily and quickly
test with varying model configurations such as extending contextual (and
non-contextual) word embeddings <cit.>
with linguistic token level
features <cit.>, and join layer document
features <cit.>
using easy to write configuration with little to no code.

What sets  apart from other frameworks is its capability of reproducing
results, efficient mini-batch creation for feature swapping for model
comparisons, and an emphasis on vectorization of natural language text
providing zero coding  construction.  The framework was written with
researchers, science related outcomers, and students in mind.

The framework's source code and installable libraries are released under the
, which is available both on GitHub and as Python  packages
along with extensive and in depth documentation, tutorials and examples.  The  documentation is fully hyperlinked, includes overview
documentation, class diagrams, and tutorials.  The framework is validated with
236 unit tests and six integration tests, most of which are automated using
continuous integration testing for both functionality and reproducibility.

prevPrevious Work

Popular  frameworks such as  have a dashboard that provides metrics,
such as training and validation loss.  However, these general purpose
frameworks offer basic performance metrics and do not provide a means of
producing higher abstraction level  specific models.  More specifically,
frameworks such as Keras, supply a very coarse  allowing solely for
cookie-cutter models.  They lack the ability to easily create and evaluate
models past this surface interface.

Frameworks such as , which are more common in academia, provide a more
straightforward simple  that is similar to the core  libraries, and
thus have the same shortcomings as a tool to bridge the gap between pure
research and reproducibility.

<cit.> is a flexible configuration
driven framework that provides ease of construction of  architectures, and thus, is the closest framework to ours.  However, it does
not have fast feature swapping persist and batch creation
capability, and lacks most of the components necessary to .

Popular packages providing support for transformer architectures such as
BERT <cit.> include .  However, this
framework only provides transformer models for contextual word embeddings.

Like the DeepDIVA <cit.>,  is
written in Python and utilizes, but does not replace, .  The goal of our
framework is:

* Reproducibility of results reproduce,

* Efficiently create and load mini-batches persist,

* Decouple the process of vectorizing data for reuse in  architectures

* Provide language specific vectorization nlfeat and  layers

All random state, including utility libraries, scientific libraries, the library, and GPU state, is consistent across each run of a Python interpreter
execution of the model's training, evaluation and testing.  Results consistency
is retained by saving this random state when saving the model, then
retrieving and resetting it after loading the model.

The order of mini-batches, and their constituent data can affect the model
performance as an aspect of training or the results of validation and testing.
This performance inconsistency is addressed by recording the order of all
data[Regardless of any user given data pre-processing or shuffling.]
and tracking the training, validation and test data splits.  Not only are
mini-batches given in the same order, the ordering in each mini-batch is also
preserved.  These dataset partitions and their order is saved to the file
system so the community has the option of distributing it along with the source
code for later experiment duplication.

In addition, the framework also saves the configuration used to recreate the
same in memory state along with the model.  This duplicates the model
structure, parameters, hyper-parameters and all other train-time memory during

stackTechnology Stack

Each “layer” of the stack builds on more general libraries to reduce the
installation footprint based on the needs of each use case.  Each library
contains the requirements for dependent third-party and lower tier packages.
The framework consists of the following libraries (see ):
The Python Library Stack

* : Utility library command line parsing, persistence and a
like inversion of
control <cit.> configuration system.

* : Parses natural language text using , generates and
vectorizes features.

* : General purpose  , much like DeepDIVA, providing
mini-batching, vectorization, and training, validation and testing of a

* : Contains language vectorization, such as word embeddings
layers, part of speech tags, named entity recognition, and head
dependencies <cit.>.

The framework provides both a command line and a  interface to
train, test and predict.  A “glue”  is used to make a  class
a dynamically generated command line with help usage message documentation.  A
set of default application classes are available with the framework, but they
can be extended to include project specific actions.  The default application
set provides interactive early stopping or epoch resetting during training.

The command line and  both use a common facade interface to the
model itself, which is conducive as an entry point to both larger projects or
simple run scripts.  However, the  interface provides evaluation
training and validation loss plots loss-snapshot.

loss-snapshotValidation and Training Loss
Plot from the NER Token Classification Application

Both interfaces provide a debugging mode that outputs a step of the model
training with batch composition, layer names, dimension calculations, using the
Python logging system, which is filterable by module or component.

The  framework provides easily configurable components to digitally
vectorize data, which in our framework, is encapsulated in a
vectorizer, which takes a particular input data and outputs a tensor.

The  library is general purpose with respect to any kind of
data, such as: text, images, audio.  Examples of data it can vectorize include
data frames, enumerated nominal values (typically used for labels),
and hot encoded vectors.

The  library provides a higher abstraction that parses natural
language text, sentence chunks, and vectorizes linguistic features.  These
vectorizers fall into one of the following categories:

token Features taken from each token, shape congruent with the number of
tokens, typically concatenated with the embedding layer include   features such as  tags,  tags, dependency tree tags and the depth
of a token in its head dependency tree.

document Features taken from the document level, typically added to a
join layer such as count sums of  parsed features.

multi-document Aggregating and shared features between
more than one document, such as overlapping  or  tags.

embedding Vectorizes text into word embeddings, such as sentence or
document text.

See nlfeat for more information on  specific feature generation.

batch-encode-tagParse Sentence Chunk and POS Tagged

For example, suppose the following sentences are to be vectorized: “The boy
hit the ball.  He did it well.”.  First the  library is used to
parse, chunk and  tag the sentence (see batch-encode-tag), then
tags are converted to one-hot encoded vectors by the  feature
vectorizer (see batch-encode-vectorize-2).

batch-encode-vectorize-2Vectorize Language

persistBatching and Persistence

Most s expect and perform well using
mini-batches <cit.> as input accept or
require batched input.  A naïve approach to generating these mini-batches
would be to re-parse and re-vectorize the data for each mini-batch over each
epoch, which is inefficient.  Many projects address this inefficiency by
pre-processing the data before training.  However, this leads to a brittle and
difficult to reproduce dataset generated set of ad-hoc text processing
scripts that are challenging re-execute, and thus, reproduce performance

A cleaner and more efficient method is to wrap this process in the framework
and create a file system scheme with the intermediate files in a configured
location so as to not clutter the project workspace, which is supported by

Another desirable feature of any framework is to easily swap in and out feature
sets and compare performance metrics, which usually takes the form of the
following steps:

* decide which features to use at train time,

* train and evaluate,

* test,

* evaluate performance,

* choose a different feature set,

* go to step a.

This incremental process highlights the need to efficiently create mini-batches.

A key observation is that each mini-batch is independent[There are
exceptions for some algorithms that need to index and fit the corpus before
vectorization.].  The  leverages the fact that mini-batches
are independent and fit nicely as independent units of work by segmenting
datasets into smaller chunks, vectorizing each chunk in parallel sub
processes, and creating one or more batches independently across each sub
processes.  This process by which data is written to the file system in a
format that is fast to reassemble is called batch encoding, and
accomplished by:

* Collecting data needed to vectorize,

* Chunking data in to equal size parts,

* Forking processes using the python multiprocessing api,

* Vectorizing each chunk in each sub-process by:

* Recreating parent memory by using configuration factory in the

* Vectorize each data chunk as separate feature sets,

* Groups the vectorized in to bundles as mini-batches, but in separate files,

* Vectorized data is almost always ready-to-go tensors.

Batch decoding is the process by which data is grouped for training as
mini-batches and is accomplished by:

* Choosing a feature set for a training run,

* Reassembling features by mini-batch, then feature as a two level directory
structure (see batch-reassembly),

* Decode each mini-batch in to a tensor, if not unserialized from the file
system as a tensor (see batch-decode),

* Copy tensors to the GPU if available,

* Cache tensors in CPU or GPU memory[Cached resources are tracked
so GPU memory is maintained.].

batch-reassemblyBatch Reassembly Process

Reassembling mini-batches by feature greatly reduces load time and memory
space, which speeds up model training and allows for more complex models.  This
leverage is most apparent when comparing pre-generated frozen large BERT model
embeddings for frozen transformers compared to a trainable model.  In the case
of the former, large data files with output tensors are read back in compared
to word piece embeddings <cit.> for a trainable

After mini-batch encoding is complete, several feature combinations can be created
in configuration, then trained, validated and tested offline.  Utility methods
exist to aggregate results in tabular form for reporting.

batch-decodeBatch Decode

nlfeatNatural Language Features

[5in]deepnlp-embeddingEmbedding Concatenation

What sets  apart from other frameworks is not only efficient mini-batch
creation and feature swapping (persist), but the tight coupling of
natural language features with a deep learning . Specifically,
vectorization of natural language features is at the heart of the utility of
the  framework, and addresses a need that is otherwise lacking in other

One such powerful capability is the concatenation of any vectorized data to
word embeddings, which is available for both contextual embeddings such as
BERT <cit.> and non-contextual embeddings such as
GLoVE <cit.> layers for supported
word embedding types.

In addition to concatenation of word embeddings, document level features can be
added to a join layer (see deepnlp-embedding).

The framework provides many layer implementations, which extend from the  class, thus any  module can be used in the framework.  To
this end, it uses layers such as the  conditional random field
implementation to create an end-to-end model for sequence classification.
Layers are configured in memory to offload the construction details to the
framework.  Other layers provided include, but are not limited to:

* BiLSTM CRF for applications such as sequence
tagging <cit.> as an end-to-end model,
which requires no coding with the exception of mapping data input to

* BERT transformer models for sentence and token classification.

* 1D convolution  that provides calculation for an arbitrarily deep
network with input and output dimensionality calculation, pooling,
mini-batch <cit.> centering and

* Expanding or contracting  feed forward linear networks with repeats
with input and output feature calculation.

* Word embedding layer capable of concatenating linguistic features at both
the token and document level vec.  Supported embeddings include
Bert <cit.>,
word2vec <cit.>,
GLoVE <cit.> and
FastText <cit.> nlfeat.

* Document level  <cit.>, and term
frequency with inverse document frequency
weighting <cit.> are also
available out of the box with no coding required

HuggingFace transformer features are available as an embeddings vectorizer and
document, sentence and token classification are available as layers.  In
addition a vectorizer is provided to map linguistic features created by
to word piece embeddings <cit.> for
concatenation with the last hidden state of the transformer.

limitsLimitations and Future Work

While the framework is thoroughly tested in the areas it was designed for, some
work remains to enable more variations in  architecture for experimentation.

Any of the  are available, but only
BERT <cit.>,
RoBERTa <cit.> and
DistilBERT <cit.> have been tested and next
sentence and masking prediction is not yet implemented.

A planned future work is to integrate the framework with TensorFlow's ,
which provides real-time graphing of metrics such as training and validation

The  framework has been presented as a viable solution to easily create
specific models with  and analysis tools to produce consistent
results.  Such frameworks are not only necessary, but vital in order to ensure
the legitimacy of the area of  in  by providing the means necessary
to produce reliable reproducible results.


# Deep learning the nonclassicality within quasi-distribution representations
from marginals

Hong-Bin Chen  Department of Engineering Science, National Cheng Kung
University, Tainan 701401, Taiwan  Center for Quantum Frontiers of Research &
Technology, NCKU, Tainan 701401, Taiwan  Physics Division, National Center for
Theoretical Sciences, Taipei 10617, Taiwan  Cheng-Hua Liu  Department of
Engineering Science, National Cheng Kung University, Tainan 701401, Taiwan
Kuan-Lun Lai  Department of Engineering Science, National Cheng Kung
University, Tainan 701401, Taiwan  Physics Division, National Center for
Theoretical Sciences, Taipei 10617, Taiwan  Bor-Yann Tseng  Department of
Engineering Science, National Cheng Kung University, Tainan 701401, Taiwan
Ping-Yuan Lo  Department of Electrophysics, National Yang Ming Chiao Tung
University, Hsinchu 300093, Taiwan  Yueh-Nan Chen  Department of Physics,
National Cheng Kung University, Tainan 701401, Taiwan  Center for Quantum
Frontiers of Research & Technology, NCKU, Tainan 701401, Taiwan  Physics
Division, National Center for Theoretical Sciences, Taipei 10617, Taiwan  Chi-
Hua Yu<EMAIL_ADDRESS>Department of Engineering Science, National
Cheng Kung University, Tainan 701401, Taiwan

###### Abstract

To unequivocally distinguish the genuine quantumness from classicality, a
widely adopted approach appeals to the negativity within a join quasi-
distribution representation as a compelling evidence for the nonclassical
essence. However, to construct a joint quasi-distribution with negativity from
experimental data typically proves to be highly cumbersome. Here we propose a
computational approach utilizing a deep generative model integrated with color
mapping to construct the bivariate joint quasi-distribution functions by
processing three marginals. We first apply our model to predict the Wigner
functions subject to thermal noises. Our model successfully predicts the
Wigner functions with a prominent accuracy by processing three marginals of
probability distributions. We also tackle a challenging problem of the
canonical Hamiltonian ensemble representation (CHER), which is developed for
characterizing the dynamical process nonclassicality. Furthermore, we also
design optimal synthetic datasets to train the model for overcoming the
ground-truth deficiency of the CHER problem. While trained with synthetic
data, the physics-informed optimization enables our model to capture the
detrimental effect of the thermal fluctuations on nonclassicality. Our
approach also provides a significant reduction of the experimental efforts of
constructing the Wigner functions of quantum states.

## Introduction

Since the birth of quantum theory, the obscure boundaries between quantum and
classical realms have always been a subject of fascination, spurring extensive
research interest Ballentine (1970); Zurek (2003); Schlosshauer (2005); Modi
_et al._ (2012). One widely adopted approach to distinctly prove genuine
quantumness as separate from classicality hinges on a classical strategy’s
inability to mimic the statistics of an experiment. The underlying philosophy
of such approach recognizes that, while a certain perspective might be
intuitively true in the classical realm, it could lead to contradictory
outcomes in quantum experiments. As such, when a classical strategy falls
short, it strongly evidences the distinct quantum essence referred to as
nonclassicality. This is especially evident in the field of quantum
information science. For instance, the renowned example of the experimental
violation Aspect  _et al._ (1981); Hensen  _et al._ (2015); Giustina  _et al._
(2015); Shalm  _et al._ (2015) of Bell’s inequality Bell (1964) accentuates
the breakdown of the EPR paradox Einstein  _et al._ (1935) formulated on the
tenets of realism and locality. This specific kind of nonclassical correlation
resulting in the violation of Bell’s inequality is termed Bell nonlocality
Brunner  _et al._ (2014).

In line with this philosophy, the inadequacy of using probability
distributions to describe quantum systems indicates the nonclassicality of
quantum states. For continuous-variable quantum systems like nano-mechanical
resonators or cavity fields, quantum states are described by quasi-
distribution functions such as the Wigner function Wigner (1932); Kenfack and
Życzkowski (2004) or the Glauber-Sudarshan $P$ representation Glauber (1963);
Sudarshan (1963). These can be interpreted as the phase-space representation
of quantum states. In the presence of strong quantum coherence, the quasi-
distribution function exhibits interference patterns, yielding negativity in
phase space. Hence, such a quantum state lacks a classical analog,
highlighting its inherent quantum essence.

This philosophy can also be employed to signify the nonclassicality of quantum
dynamical processes Lambert  _et al._ (2010); Rahimi-Keshari  _et al._ (2013);
O’Reilly and Olaya-Castro (2014); Sabapathy (2016); Smirne  _et al._ (2019);
Milz  _et al._ (2020); Seif  _et al._ (2022) through the newly emerging
canonical Hamiltonian ensemble representation (CHER) Chen  _et al._ (2019);
Chen and Chen (2021); Lin  _et al._ (2022). In CHER, the dynamical behaviors
of open quantum systems are described by the ensemble average of unitary
evolutions characterized by a Hamiltonian ensemble (HE) Kropf  _et al._
(2016); Chen  _et al._ (2018); Chen (2022), which comprises a collection of
parameterized system Hamiltonian operators and the corresponding quasi-
distribution function. Crucially, the presence of negativity in the quasi-
distribution function indicates the establishment of system-environment
quantum correlations, thereby emphasizing the quantum essence of the dynamical
process. Consequently, the quasi-distribution representation emerges as a
cornerstone in characterizing the inherent nonclassicality in quantum
dynamical processes.

Despite its potential in characterizing the nonclassicality, constructing a
multivariate joint quasi-distribution with negativity from experimental data
remains a cumbersome task. For instance, to construct the Wigner function, one
might first determine the complete quantum state through methods such as
tomographically measuring rotated quadratures Smithey  _et al._ (1993);
Breitenbach  _et al._ (1997); Lvovsky  _et al._ (2001), transformations from
different quasi-distribution representations Deléglise  _et al._ (2008), or
characteristic-function tomography Flühmann and Home (2020). Afterward, the
Wigner function can be constructed by post-processing the measured data.
Alternatively, one can also directly scan the Wigner function point-by-point
in phase space Lutterbach and Davidovich (1997); Banaszek  _et al._ (1999);
Bertet  _et al._ (2002); Vlastakis  _et al._ (2013); Winkelmann  _et al._
(2022). All these approaches demand extensive experimental efforts. Even in
the machine-learning-assisted method Ahmed  _et al._ (2021), thousands of data
points are still required.

The complexity escalates when dealing with the nonclassicality of dynamical
processes, as it involves quantum process tomography. According to the
underlying Lie algebraic structure, only specific marginals can be efficiently
solved, suggesting that the CHER can only be built from several marginal
distributions. This is reminiscent of a very challenging mathematical issue
that has been extensively studied Cuadras  _et al._ (2002); Lin  _et al._
(2014); Takeuchi (2010); Cicalese  _et al._ (2019); Ong  _et al._ (2021).
Particularly, in the context of CHER, complications arise as the quasi-
distribution may exhibit negative values, rendering previously developed
mathematical techniques unreliable. Therefore, a reliable novel technique to
construct quasi-distribution functions, allowing for negativity, from given
marginal distributions is highly demanded. Such a technique would not only be
crucial for constructing CHER, but might also significantly cut down the
experimental efforts needed for constructing the Wigner functions of quantum
states.

Here we propose to harness the power of deep generative models (DGMs) to
generate the corresponding joint quasi-distribution by processing a series of
marginals (Fig. 1). DGMs excel in discerning hidden patterns in expansive
datasets, proving invaluable in extracting new insights. Our DGM integrates
the ResNet structure He  _et al._ (2016) to seamlessly generate the bivariate
joint quasi-distribution for both challenges. To further enhance our DGM’s
pattern detection within an image, we have devised a color mapping strategy to
translate the joint quasi-distribution into three monochromatic images. To
validate our approach and to offer an intuitive result presentation, we
present examples including the Wigner function of quantum harmonic oscillator
exposed to thermal bath and the CHER of qubit-pair pure dephasing at various
temperatures.

Our appraoch offers a streamlined computational alternative, lessening the
rigorous experimental demands of Wigner function measurements. Once the well-
trained DGM receives three marginals, which correspond to three probability
distributions in real or momentum space, it can generate the desired Wigner
function. Simultaneously, we have addressed the issue of lacking ground truth
(GT) for predicting the CHER. By leaning on prior knowledge about the target
marginals, synthetic training datasets can be effectively optimized. We
demonstrate that our physics-informed approach adeptly constructs high-
precision Wigner functions and identifies the detrimental effect of thermal
fluctuation on the nonclassicality, especially as temperature rises in the
qubit-pair pure dephasing dynamics.

### Constructing joint quasi-distribution functions from marginals

Figure 1: Construction of bivariate joint quasi-distribution from a set of
marginals. While the joint quasi-distribution function offers invaluable
insights into the quantum essence, its construction generically requires
substantial efforts. In contrast, constructing marginals is relatively
straightforward. Here we develop a DGM integrating color mapping. This DGM
enables the construction of bivariate joint quasi-distribution functions by
processing three marginals. These marginals can originates from experimental
raw data or be theoretically derived from a prototypical physical model.

The quasi-distribution representation is a prevalent approach in quantum
physics. Its negativity conveys a visualized insight into genuine quantum
characteristics. However, in many practical scenarios, its construction can be
challenging or demands substantial efforts. For instance, the direct point-by-
point measurement of Wigner functions in phase space is experimentally
intensive, and the construction of the CHERs for higher-dimensional dynamical
processes is typically infeasible. In contrast, the marginals derived from the
joint quasi-distribution functions describing nonclassicality are generically
standard probability distributions. This renders the marginals also physically
meaningful and, most importantly, makes them accessible for construction using
experimental raw data.

Motivated by this observation, our goal is to devise an approach for
constructing joint quasi-distribution functions, which permit negativity,
using three marginals. As illustrated in Fig. 1, we harness artificial
intelligence techniques, specifically designing a DGM integrated with color
mapping. Once we have the marginals, whether sourced from experimental data or
theoretical calculations, our well-trained DGM efficiently constructs the
corresponding joint quasi-distribution function with notable precision.

### Wigner functions and the marginals

The motional state of a continuous-variable quantum system can be equivalently
described by either the density matrix $\rho$ or the Wigner function
$\mathcal{W}(x,p)=(\pi\hbar)^{-1}\int{\langle
x+x^{\prime}|}\rho{|x-x^{\prime}\rangle}\exp(-i2px^{\prime}/\hbar)dx^{\prime}$,
which is a real-valued quasi-distribution function permitting negativity. Such
a phase-space representation has been proven useful for visualizing quantum
states. Particularly, the negativity indicates the inherent quantum essence of
the states. Most importantly, a unique property of the Wigner function is that
the marginals are standard probability distributions over certain quadratures,
and therefore easily experimentally measurable, in contrast to the point-by-
point measurement of the Wigner function in phase space Lutterbach and
Davidovich (1997); Banaszek  _et al._ (1999); Bertet  _et al._ (2002);
Vlastakis  _et al._ (2013); Winkelmann  _et al._ (2022). However, to construct
a joint quasi-distribution from a few number of specific marginals remains a
challenging mathematical issue Cuadras  _et al._ (2002); Lin  _et al._ (2014);
Takeuchi (2010). Moreover, the well-developed algorithm of inverse Radon
transform requires an extensive measurement over rotated quadratures Smithey
_et al._ (1993); Breitenbach  _et al._ (1997); Lvovsky  _et al._ (2001) for
tomographically constructing the Wigner function.

We design a DGM for constructing the Wigner function using merely three
marginals. To showcase our approach, here we specifically consider a quantum
harmonic oscillator in coherent or cat sates exposed to thermal bath. The two
natural marginals, $\mathcal{W}(x)$ and $\mathcal{W}(p)$, possess clear
physical interpretation as the probability distributions in real and momentum
space. While the third one $\mathcal{W}(u)$ over the oblique variable
$u=(x+p)/\sqrt{2}$ is accessible in real space after a $\pi/4$-rotation of the
quantum states (Supplementary Note 1).

### Data representation and color mapping

Figure 2: Generation of training data with color mapping. a In the case of the
Wigner function, the training dataset can be generated efficiently. Here we
exemplify our approach by specifically generating parameterized coherent
states and cat states subject to thermal noises. b, c The color mapping
consists of three functions mapping the (rescaled) height $\zeta$ of the
generated training data to a set of RGB values, leading to three monochromatic
images to be learned by the DGM; meanwhile, it can also be merged into a
colored image for visualizing the data. d The three marginals
$(\mathcal{W}(x),\mathcal{W}(p),\mathcal{W}(u))$ can be efficiently derived
from the Wigner function $\mathcal{W}(x,p)$ by integral over redundant
variable. e Each single marginal is sampled into 721 pixels and each
monochromatic image is sampled into a $256\times 256$ pixel array. Finally,
the three marginals are combined with the three monochromatic images,
constituting the feature and the label of a single training datum,
respectively. Here we show the raw data structure fed into our deep generative
model and the corresponding visualization.

To train a supervised DGM, it is necessary to prepare a labelled training
dataset. In the case of the Wigner functions, the training data can be
generated efficiently (Methods). The data representation is schematically
illustrated in Fig. 2. The feature of each training datum consists of three
marginals $(\mathcal{W}(x),\mathcal{W}(p),\mathcal{W}(u))$, which can be
derived analytically (Fig. 2d), even in the presence of thermal noises. Each
marginal is numerically sampled into 721 pixels as the raw data. To take the
advantage of our DGM in discerning patterns and relationships within the RGB
channels of an image, each generated Wigner function $z=\mathcal{W}(x,p)$ is
translated into three monochromatic images according to a color mapping (Figs.
2b and c). The color mapping comprises three functions
$(f_{\mathrm{R}}(\zeta),f_{\mathrm{G}}(\zeta),f_{\mathrm{B}}(\zeta))$, mapping
the (rescaled) height value $\zeta=(z+0.45)/0.9$ to a set of RGB values,
resulting in the three monochromatic images. Note that the color mapping
should be optimized according to the training data to be learned (Methods).
Finally, the output monochromatic images serve as the label to be learned by
the DGM, each of which is sampled into a $256\times 256$ pixel array. The
visualization of a single training datum is shown in Fig. 2e (Supplementary
Note 3). Crucially, the color mapping associates the positivity and the
negativity in the generated data to the R- and B-channels, respectively. This
helps our DGM pay particular attention to the negativity in the Wigner
function, capturing the nonclassicality of quantum states.

### ResNet-based deep generative model

Our DGM is structured into six stages, each constructed by repeatedly stacking
two core building blocks: the identity and deconvolution blocks. As our model
aims to produce three monochromatic images from input marginals, the main
stream of both blocks utilizes three deconvolutional layers. These layers
extract key patterns from the input marginals and expand them to produce the
output monochromatic images of a quasi-distribution. For the model to handle
this intricate task, it requires more depth and additional deconvolutional
layers. However, this intensifies the vanishing gradient problem, which
significantly impedes our model’s training.

To address this challenge, we have incorporated a shortcut to each block,
known as the ResNet structure He  _et al._ (2016). Thanks to the skip
connections in the ResNet structure, information from previous layers is
retained, substantially mitigating the vanishing gradient problem. This allows
us to deepen our model, enhancing its capabilities. Further insights and
comprehensive layout of our DGM are shown in Supplementary Notes 4.

### Constructed Wigner functions by the DGM

Figure 3: Verifying the DGM’s performance on the Wigner functions. The top row
shows several GT Wigner functions sampled from the testing data used to derive
the marginals for the well-trained DGM. The middle row shows the corresponding
predictions. A close comparison reveals the DGM’s impressive accuracy, with
differences being almost indistinguishable. The bottom row depicts the pixel-
wise absolute differences in grayscale. The DGM precisely predicts noisy
coherent states with notable precision. Prominent numerical discrepancies only
emerges in regions showcasing significant interference patterns in the noisy
cat states.

After training our DGM, we apply it to predict the Wigner functions from the
three testing data marginals. Figure 3 shows several representative samples.
The top row shows the GT Wigner functions used to derive the marginals input
into the DGM. The middle row shows the DGM’s predictions. A close comparison
underscores the DGM’s accuracy with negligible differences between GT and
predictions. The bottom row highlights the pixel-wise absolute difference
between GT and predictions in grayscale. The DGM’s predictions are notably
precise for noisy coherent states, as seen in the first and fifth columns.
However, regions with significant interference patterns in the noisy cat
states show more pronounced differences. Further comprehensive quantitative
assessment of the model’s performance is presented in Supplementary Note 5.

### Dynamical process nonclassicality characterized by the CHER

To further demonstrate the potential of our DGM, we apply our model to the
problem of CHER, which is even more challenging due to its GT-deficiency.

In the theory of dynamical process nonclassicality, the mathematical toolbox
of the HE Kropf  _et al._ (2016); Chen (2022)
$\\{(p_{\lambda},\widehat{H}_{\lambda})\\}_{\lambda}$ is leveraged to
formulate the classical strategy simulating an incoherent dynamics. A HE is a
collection of traceless Hermitian operators
$\widehat{H}_{\lambda}\in\mathfrak{su}(n)$ associated with a probability
$p_{\lambda}$ of occurrence. It has been shown that if a system can merely
establish classical correlations with its environment during their
interaction, then the resulting reduced dynamics $\mathcal{E}_{t}$ of the
system admits a HE decomposition Chen  _et al._ (2018) written as

$\mathcal{E}_{t}\\{\rho(0)\\}=\int
p_{\lambda}\widehat{U}_{\lambda}(t)\rho(0)\widehat{U}^{{\dagger}}_{\lambda}(t)d\lambda$
(1)

for a certain HE, $\widehat{U}_{\lambda}(t)=\exp(-i\widehat{H}_{\lambda}t)$,
and the initial state $\rho(0)$ of the system. Therefore, irrespective of the
inaccessible environmental degrees of freedom in the vast majority of
authentic experimental setups, we classify an incoherent dynamics
$\mathcal{E}_{t}$ as classical if it admits a HE-simulation (1). An intuition
of this criterion can be acquired from the observation that Eq. (1) can be
explained classically as a statistical mixture of random unitary channels. On
the contrary, $\mathcal{E}_{t}$ is nonclassical if it admits no HE-
simulations. Such nonclassical dynamics arises from the consumption of
nonclassical correlations, not reproducible by classical resources. The
nonexistence of simulating HEs can be proven by the necessity to resort to a
nonclassical HE encapsulating a quasi-distribution $\wp_{\lambda}$ showing
negative values. Consequently the negativity in $\wp_{\lambda}$ is an
indicator of the nonclassicality of $\mathcal{E}_{t}$.

To further promote $\wp_{\lambda}$ to a characteristic representation over the
frequency domain for $\mathcal{E}_{t}$, referred to as CHER, the ensemble-
averaged simulation (1) has been recast into the FToG formalism Chen  _et al._
(2019); Chen and Chen (2021)

$\mathcal{E}^{(\widetilde{L})}_{t}=\int_{\mathcal{G}}\wp_{\lambda}e^{-i\lambda\widetilde{L}t}d\lambda,$
(2)

where $\mathcal{E}^{(\widetilde{L})}_{t}$ is the linear map of
$\mathcal{E}_{t}$ expressed with respect to the adjoint representation
$\\{\widetilde{L}\\}$ of the generators of $\mathfrak{u}(n)$. Noteworthily,
this formalism is derived according to the Lie-algebraic structure underlying
the HE, rather than the conventional Fourier transform. The FToG (2)
associates the quasi-distribution $\wp_{\lambda}$ with the dynamical process
$\mathcal{E}_{t}$, i.e.,
$\wp_{\lambda}\mapsto\mathcal{E}_{t}^{(\widetilde{L})}$, highlighting the role
of $\wp_{\lambda}$ as a CHER for $\mathcal{E}_{t}$. More specifically, to
explore the HE-simulation for $\mathcal{E}_{t}$ means to construct the CHER
$\wp_{\lambda}$ of a joint quasi-distribution by solving the FToG (2). For a
given pure dephasing dynamics $\mathcal{E}_{t}$ of any dimension, the Lie-
algebraic structure of the FToG leads to the marginals along the root vectors
of the Lie algebra (Supplementary Note 7). Then the corresponding CHER
$\wp_{\lambda}$ can only be built from the marginals solved from the FToG,
rendering the direct construction of the CHER intractable.

To compensate this deficiency, we apply the DGM to predict the CHER by
processing a set of marginals solved from the FToG. To facilitate an intuitive
result presentation, here we specifically consider a two-dimensional problem
over the $x_{1}$-$x_{13}$ plane, as well as the $u=(x_{1}+x_{13})/\sqrt{2}$
coordinate describing the correlation between them. We stress that our
approach can be generalized to higher dimension by appropriately including
more marginals and correlations.

### Preparation of synthetic training data

Unlike the situation with the Wigner function, the challenge of CHER stems
from its GT-deficiency. Specifically, the absence of explicit CHER solution
from the FToG (2) makes the generation of adequate training data problematic.
To address this issue, we opt to generate synthetic training datasets. In
these datasets, the joint quasi-distribution is crafted from a combination of
one bivariate Gaussians with a pair of Gaussians with opposite peaks,
simulating potential negativity in the CHER. This allows us to efficiently
derive the marginals. Additionally, the synthetic datasets are informed by
physical principles, making them adaptable. By tweaking relevant parameters,
they can be optimized for sufficiently reflecting the target marginals of a
specific quantum dynamics model (Methods and Supplementary Note 8). This
optimization aids in enhancing the proficiency and precision of our DGM. Once
the synthetic data are generated, they can be encoded into the raw data
following the process illustrated in Fig. 2.

Furthermore, we set our sights on training the DGM to forecast two separate
quantum dynamics models. To achieve this, we generate two unique synthetic
training datasets, each optimized for specific quantum dynamics models. This
optimization process begins by examining the profile of the marginals solved
from the FToG (2) for the intended quantum models. Armed with this knowledge,
we can fine-tune the parameters during the synthetic data generation. This
ensures that our synthetic marginals align closely with the nuances of the
quantum model to be solved (Supplementary Note 8).

### Quantum pure dephasing dynamics

As a concrete paradigm, we consider an extended spin-boson model consisting of
a non-interacting qubit pair coupled to a common boson bath with total
Hamiltonian
$\widehat{H}_{\mathrm{T}}=\sum_{j=1,2}\omega_{j}\hat{\sigma}_{z,j}/2+\sum_{\mathbf{k}}\omega_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}+\sum_{j,\mathbf{k}}\hat{\sigma}_{z,j}\otimes(g_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{\dagger}+g_{\mathbf{k}}^{\ast}\hat{b}_{\mathbf{k}})$.
The qubit pair undergoes a pure dephasing dynamics characterized by four
dephasing factors and the corresponding CHER
$\wp(x_{1},x_{6},x_{13})=\wp_{1,13}(x_{1},x_{13})\wp_{6}(x_{6})$ is governed
by the FToG (2) explicitly written as:

$\left\\{\begin{array}[]{l}\phi_{1}(t)=\exp[i\vartheta(t)-\Phi(t)]=\int_{\mathbb{R}}\wp_{1}(x_{1})\exp[-ix_{1}t]dx_{1}\\\
\phi_{9}(t)=\exp[-4\Phi(t)]=\int_{\mathbb{R}^{2}}\wp_{1,13}(x_{1},x_{13})\exp[-ix_{1}t]\exp[-ix_{13}t]dx_{1}dx_{13}\\\
\phi_{13}(t)=\exp[-i\vartheta(t)-\Phi(t)]=\int_{\mathbb{R}}\wp_{13}(x_{13})\exp[-ix_{13}t]dx_{13}\\\
\phi_{6}(t)=1=\int_{\mathbb{R}}\wp_{6}(x_{6})\exp[-ix_{6}t]dx_{6}\end{array}\right.,$
(3)

with $\vartheta(t)=4\int_{0}^{\infty}[\mathcal{J}(\omega)/\omega^{2}](\omega
t-\sin\omega t)d\omega$,
$\Phi(t)=4\int_{0}^{\infty}[\mathcal{J}(\omega)/\omega^{2}]\coth(\hbar\omega/2k_{\mathrm{B}}T)(1-\cos\omega
t)d\omega$, and
$\mathcal{J}(\omega)=\sum_{\vec{k}}|g_{\vec{k}}|^{2}\delta(\omega-\omega_{\vec{k}})$
being the spectral density. Note that the first three lines determine the
three marginals along the $x_{1}$-, $u$-, and $x_{13}$-axis, respectively,
which will be fed into a well-trained model to forecast a
$\wp_{1,13}(x_{1},x_{13})$; while the last line leads to an independent
component $\wp_{6}(x_{6})=\delta(x_{6})$. In the following, we demonstrate two
types of spectral density, i.e., the super-Ohmic family and the Drude-Lorentz
spectral densities, at various temperatures. Further details are shown in
Supplementary Notes 9.

### Verifying model performance

After training with our synthetic data, we proceed to verify the model’s
performance on this data, using a verification procedure consistent with the
one used for the Wigner functions. We also compare both the joint quasi-
distributions and the derived marginals of the GT and the model’s prediction.
This allows us to affirm the reliability of our well-trained models on both
synthetic datasets. Further details are shown in Supplementary Notes 10.

Nevertheless, when forecasting the CHER of a dynamics, the GT-deficiency makes
the direct comparison of the joint quasi-distributions infeasible.
Subsequently, we show the predicted CHERs for both quantum spectral densities
and estimate the pixel-averaged $\mathrm{L}_{1}$ norm, defined as
$\sum_{\mathrm{pixel}}|\mathrm{GT}-\mathrm{Prediction}|/721$, for marginals
derived both from FToG and the predicted CHERs. This methodology allows us to
gauge the reliability of our model in forecasting the CHER, even with the
obstacles of GT-deficiency.

### Prediction of CHER and nonclassicality

Figure 4: CHERs predicted by the DGM for the super-Ohmic spectral density. The
first row shows the predicted CHER for the super-Ohmic spectral density
(Ohmicity $s=2$) with increasing temperature $T$. The central peak is
suppressed by the increasing thermal fluctuations in the boson bath. The
shallow negative regions are magnified in the second row for clarity,
indicating the nonclassicality of the dynamics. The negative region is also
gradually attenuated with increasing temperature, reflecting the detrimental
effects of thermal fluctuation on the nonclassical essence. The comparison
between the marginals derived from the predicted CHERs (red solid curves) and
those of GT (dashed black curves) are shown in the last three rows with
estimated pixel-averaged $\mathrm{L}_{1}$ norm shown at the upper-right corner
of each panel. These comparisons affirm the reliability of the predicted
CHERs. Figure 5: CHERs predicted by the DGM for the Drude-Lorentz spectral
density. The first row shows the predicted CHER for the Drude-Lorentz spectral
density with increasing temperature $T$. Similar suppression of the central
peak with increasing $T$ can also be observed. The shallow negative regions
are magnified in the second row for clarity. The phenomenon of quantum-to-
classical transition is revealed by the negative region vanishing with
increasing $T$. The comparison between the marginals derived from the
predicted CHERs (red solid curves) and those of GT (dashed black curves) are
shown in the last three rows with estimated pixel-averaged $\mathrm{L}_{1}$
norm shown at the upper-right corner of each panel, which underpin the
reliability of our model.

Figure 4 shows the prediction for the quantum dynamics models (3) with a
specific super-Ohmic spectral density (Ohmicity $s=2$) in the family. The
first row shows the predicted CHER with increasing temperature $T$, which
suppresses the central peak of the CHER. Each CHER possesses a shallow
negative region indicating the nonclassicality of the dynamics, which is
magnified in the second row for clarity. Crucially, the negative region is
gradually attenuated with increasing temperature, reflecting the effects of
increasing thermal fluctuation. This is in line with the intuition that the
thermal fluctuation is detrimental to the nonclassical essence of a quantum
system. Due to the lack of the GT CHER constructed from the FToG (2), to
underscore the reliability of the predictions, we have also numerically
derived the marginals (red solid curves) from the predicted CHERs and compared
with those of GT (dashed black curves) solved from the quantum dynamics models
(3), as shown in the last three rows. The pixel-averaged $\mathrm{L}_{1}$ norm
is shown at the upper-right corner of each panel. Generically, the prediction
fits the GT better at higher temperatures, resulting in lower errors.
Furthermore, the errors are two-order of magnitude lower than the peaking
values of the marginals. This affirms the reliability of the predicted CHERs.

Figure 5 shows the prediction for the Drude-Lorentz spectral density. We show
the predicted CHER in the first row with increasing temperature $T$. Similar
suppression of the central peak with increasing $T$ can also be observed. The
negative regions are magnified in the second row. The nonclassicality is
destroyed at a high enough $T$, revealing a quantum-to-classical transition
due to the thermal fluctuations. We also show the comparison between the
derived the marginals (red solid curves) from the prediction and the GT
marginals (dashed black curves) in the last three rows with estimated pixel-
averaged $\mathrm{L}_{1}$ norm shown at the upper-right corner, which underpin
the reliability of our model.

## Discussion

In this study, we’ve reached a notable milestone in quantum information
science by creatively employing deep generative models (DGMs) to craft
multivariate joint quasi-distribution functions. This is a critical step in
detailing the nonclassical nature of quantum states and their dynamics. Our
method skillfully sidesteps the constraints of traditional techniques, such as
the Wigner function, by utilizing sophisticated neural networks like the
ResNet model, along with a novel color mapping technique. This advanced
approach has not only streamlined the analysis of experimental data but has
also introduced a new computational framework in quantum physics. Our DGMs’
accuracy in predicting quantum states, even amidst thermal noise, represents a
substantial advancement in our comprehension and measurement of
nonclassicality. This work signifies a major stride in merging quantum
mechanics with artificial intelligence, setting the stage for revolutionary
progress in quantum computing and information processing.

The deployment of DGMs in our research signifies a pivotal transformation in
quantum physics, particularly in characterizing and analyzing quantum states.
These models’ superior computational abilities have surpassed conventional
boundaries, allowing for a more intricate and precise depiction of quantum
states. The capacity of the DGMs to analyze and decode complex datasets, such
as those from quantum experiments, has unveiled new pathways for grasping
quantum phenomena. By accurately forecasting the behavior of quantum systems,
inclusive of the nuances introduced by thermal noise, these models offer an
expanded comprehension of quantum nonclassicality. This breakthrough not only
boosts the precision of quantum state analysis but also considerably lowers
the experimental and computational efforts typically needed in this domain.
The integration of DGMs signals a shift in both technological progress and the
conceptual approach to quantum state analysis and representation, potentially
sparking novel discoveries and applications in quantum computing, simulation,
and data handling.

Despite the advances made in this research, there are intrinsic limitations to
our method. A key challenge is the dependence on the quality and diversity of
the training data for the DGMs. Inaccuracies or biases within this data could
result in less reliable forecasts, especially in scenarios vastly differing
from the training set. Moreover, the computational complexity of DGMs, when
handling high-dimensional quantum systems, presents a significant hurdle,
potentially affecting the scalability of our method. Another concern is the
interpretability of the model’s output; like many deep learning models, DGMs
often operate as ’black boxes’, which can obscure the rationale behind their
predictions.

To overcome these limitations, future studies should concentrate on refining
data collection and preprocessing methods, to ensure a thorough and impartial
dataset for training the DGMs. Advancing data augmentation techniques may also
contribute to developing a robust model that can generalize effectively to
new, uncharted scenarios. Addressing the computational hurdles, optimizing the
neural network structure, and employing more efficient training algorithms may
offer solutions. Moreover, the adoption of explainable AI methods could
enhance the DGMs’ interpretability, enabling researchers to delve deeper into
the model’s decision-making process. Continued collaboration between quantum
physicists and AI experts is essential to refine these models, ensuring they
stay relevant and effective in the rapidly evolving quantum research
landscape.

Looking forward, the scope for further advancements in this field is immense.
Our research illustrates how the convergence of quantum physics and artificial
intelligence can lead to innovative approaches in quantum computing,
simulation, and information handling. Future research should aim to improve
the training and scalability of DGMs to manage more intricate quantum systems,
and possibly explore new quantum phenomena. Furthermore, improving the models’
interpretability is crucial for deeper insight into their predictions and for
building trust in their applications. As quantum technology advances, the
incorporation of sophisticated computational models like DGMs will be vital in
unlocking fresh capabilities and fostering innovation. Interdisciplinary
collaboration, especially between quantum physicists and AI experts, will be
central to fully leveraging these technologies, leading to landmark
discoveries and applications in quantum science.

Our study signifies an important advance in quantum information science
through the successful integration of DGMs into the analysis and
characterization of quantum states and dynamics. This novel approach surmounts
several limitations of conventional methods, offering a more nuanced and
comprehensive understanding of quantum nonclassicality. The precision and
computational efficiency brought by our DGMs not only improve the accuracy of
quantum state characterization but also open new avenues for exploring
intricate quantum phenomena. The effective implementation of these models
underscores the vast potential of artificial intelligence in enhancing our
understanding of quantum mechanics and its practical applications.

## Methods

### Generation of Wigner function training data of noisy states

In our Wigner function training dataset, we specifically consider two types of
parameterized states, i.e., the noisy coherent states and the noisy cat
states. Since the coherent parameter $\alpha\in\mathbb{C}$ is a complex
number, its real $\mathrm{Re}[\alpha]$ and imaginary parts
$\mathrm{Im}[\alpha]$ are independently sampled in the range $[-2,2]$. While
the relative phase $\theta$ in the noisy cat state is sampled in the range
$[0,2\pi)$. The decoherence effect is characterized by $\mu\in[0.5,1]$, while
the thermal noise model Xiong  _et al._ (2015) is characterized by
$\nu=(1-\mu^{2})\bar{n}$ with $\bar{n}\in[0,2]$ denoting the averaged
excitation in the thermal bath. We have generated 16,000 noisy coherent states
and 18,000 noisy cat states in the Wigner function training dataset. Since,
for both states, the Wigner function $\mathcal{W}(x,p)$ and the three
marginals $(\mathcal{W}(x),\mathcal{W}(p),\mathcal{W}(u))$ can be derived
analytically, the Wigner function training dataset can be generated
efficiently. The detialed analytic expressions are presented in Supplementary
Note 1.

### Color mapping

To deal with the joint quasi-distribution in the sense of an image, we design
a color mapping to convert each datum into three monochromatic images. The
color mapping consists of three functions

$\left\\{\begin{array}[]{l}f_{\mathrm{R}}(\zeta)=2\times\frac{1.148}{(e^{-25(\zeta-(\zeta_{0}-0.12))}+1)(e^{5(\zeta-(\zeta_{0}+0.45))}+1)}-1\\\
f_{\mathrm{G}}(\zeta)=2\times e^{\frac{-(\zeta-\zeta_{0})^{2}}{0.0392}}-1\\\
f_{\mathrm{B}}(\zeta)=2\times\frac{1.148}{(e^{25(\zeta-(\zeta_{0}+0.12))}+1)(e^{-5(\zeta-(\zeta_{0}-0.45))}+1)}-1\end{array}\right.,$
(4)

where $\zeta_{0}$ controls the peaking positions of the RGB curves, and the
hight value of the generated data $z\mapsto\zeta$ is rescaled such that
$\zeta\in[0,1]$ for most of the training data. These parameters should be
determined according to the statistics of the training data to be learned.
Over-compressing will smear the details in the joint quasi-distributions, and
under-compressing will truncate the peaking values of the quasi-distributions.
In the case of Wigner function, we set $\zeta_{0}=1/2$ and
$\zeta=(z+0.45)/0.9$. While for the case of CHER, we set $\zeta_{0}=1/5.5$ and
$\zeta=(z+0.01)/0.055$. Further advantages of the color mapping on improving
our model and the raw data structure are discussed in Supplementary Note 3.

### Generation of synthetic training data

To train the DGM for solving the GT-deficient problem, we need to generate a
sufficient amount of relevant training data. However, it is infeasible for us
to achieve this by solving the FToG of a quantum dynamics. We circumvent this
obstacle by generating synthetic training datasets. We combine a conventional
bivariate Gaussian with a pair of Gaussians with opposite peaks according to

$\wp(x_{1},x_{13})=p(x_{1},x_{13})+Ap^{\prime}(x_{1},x_{13})-Ap^{\prime\prime}(x_{1},x_{13}),$
(5)

where $p(x_{1},x_{13})$, $p^{\prime}(x_{1},x_{13})$, and
$p^{\prime\prime}(x_{1},x_{13})$ are three randomly generated (conventional)
Gaussians with individual statistical parameters, i.e., mean value, standard
deviation, and correlation, and $A$ is a random amplitude. It is easy to see
that $\wp(x_{1},x_{13})$ is normalized. Further details are shown in
Supplementary Note 8.

## Data availability

The data that support the findings of this study are available upon reasonable
request from the corresponding authors.

## References

  * Ballentine (1970) L. E. Ballentine, “The statistical interpretation of quantum mechanics,” Rev. Mod. Phys. 42, 358 (1970).
  * Zurek (2003) W. H. Zurek, “Decoherence, einselection, and the quantum origins of the classical,” Rev. Mod. Phys. 75, 715 (2003).
  * Schlosshauer (2005) M. Schlosshauer, “Decoherence, the measurement problem, and interpretations of quantum mechanics,” Rev. Mod. Phys. 76, 1267 (2005).
  * Modi  _et al._ (2012) K. Modi, A. Brodutch, H. Cable, T. Paterek, and V. Vedral, “The classical-quantum boundary for correlations: Discord and related measures,” Rev. Mod. Phys. 84, 1655 (2012).
  * Aspect  _et al._ (1981) A. Aspect, P. Grangier, and G. Roger, “Experimental tests of realistic local theories via bell’s theorem,” Phys. Rev. Lett. 47, 460 (1981).
  * Hensen  _et al._ (2015) B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, “Loophole-free bell inequality violation using electron spins separated by 1.3 kilometres,” Nature 526, 682 (2015).
  * Giustina  _et al._ (2015) M. Giustina, M. A. M. Versteegh, S. Wengerowsky, J. Handsteiner, A. Hochrainer, K. Phelan, F. Steinlechner, J. Kofler, J.-A. Larsson, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, J. Beyer, T. Gerrits, A. E. Lita, L. K. Shalm, S. W. Nam, T. Scheidl, R. Ursin, B. Wittmann, and A. Zeilinger, “Significant-loophole-free test of bell’s theorem with entangled photons,” Phys. Rev. Lett. 115, 250401 (2015).
  * Shalm  _et al._ (2015) L. K. Shalm, E. Meyer-Scott, B. G. Christensen, P. Bierhorst, M. A. Wayne, M. J. Stevens, T. Gerrits, S. Glancy, D. R. Hamel, M. S. Allman, K. J. Coakley, S. D. Dyer, C. Hodge, A. E. Lita, V. B. Verma, C. Lambrocco, E. Tortorici, A. L. Migdall, Y. Zhang, D. R. Kumor, W. H. Farr, F. Marsili, M. D. Shaw, J. A. Stern, C. Abellán, W. Amaya, V. Pruneri, T. Jennewein, M. W. Mitchell, P. G. Kwiat, J. C. Bienfang, R. P. Mirin, E. Knill, and S. W. Nam, “Strong loophole-free test of local realism,” Phys. Rev. Lett. 115, 250402 (2015).
  * Bell (1964) J. S. Bell, “On the einstein podolsky rosen paradox,” Physics 1, 195 (1964).
  * Einstein  _et al._ (1935) A. Einstein, B. Podolsky, and N. Rosen, “Can quantum-mechanical description of physical reality be considered complete?” Phys. Rev. 47, 777 (1935).
  * Brunner  _et al._ (2014) N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, “Bell nonlocality,” Rev. Mod. Phys. 86, 419 (2014).
  * Wigner (1932) E. P. Wigner, “On the quantum correction for thermodynamic equilibrium,” Phys. Rev. 40, 749 (1932).
  * Kenfack and Życzkowski (2004) A. Kenfack and K. Życzkowski, “Negativity of the wigner function as an indicator of non-classicality,” J. Opt. B Quantum Semicl. Opt. 6, 396 (2004).
  * Glauber (1963) R. J. Glauber, “Coherent and incoherent states of the radiation field,” Phys. Rev. 131, 2766 (1963).
  * Sudarshan (1963) E. C. G. Sudarshan, “Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams,” Phys. Rev. Lett. 10, 277 (1963).
  * Lambert  _et al._ (2010) N. Lambert, C. Emary, Y.-N. Chen, and F. Nori, “Distinguishing quantum and classical transport through nanostructures,” Phys. Rev. Lett. 105, 176801 (2010).
  * Rahimi-Keshari  _et al._ (2013) S. Rahimi-Keshari, T. Kiesel, W. Vogel, S. Grandi, A. Zavatta, and M. Bellini, “Quantum process nonclassicality,” Phys. Rev. Lett. 110, 160401 (2013).
  * O’Reilly and Olaya-Castro (2014) E. J. O’Reilly and A. Olaya-Castro, “Non-classicality of the molecular vibrations assisting exciton energy transfer at room temperature,” Nat. Commun. 5, 3012 (2014).
  * Sabapathy (2016) K. K. Sabapathy, “Process output nonclassicality and nonclassicality depth of quantum-optical channels,” Phys. Rev. A 93, 042103 (2016).
  * Smirne  _et al._ (2019) A. Smirne, D. Egloff, M. G. Díaz, M. B. Plenio, and S. F. Huelga, “Coherence and non-classicality of quantum markov processes,” Quantum Sci. Technol. 4, 01LT01 (2019).
  * Milz  _et al._ (2020) S. Milz, D. Egloff, P. Taranto, T. Theurer, M. B. Plenio, A. Smirne, and S. F. Huelga, “When is a non-markovian quantum process classical?” Phys. Rev. X 10, 041049 (2020).
  * Seif  _et al._ (2022) A. Seif, Y.-X. Wang, and A. A. Clerk, “Distinguishing between quantum and classical markovian dephasing dissipation,” Phys. Rev. Lett. 128, 070402 (2022).
  * Chen  _et al._ (2019) H.-B. Chen, P.-Y. Lo, C. Gneiting, J. Bae, Y.-N. Chen, and F. Nori, “Quantifying the nonclassicality of pure dephasing,” Nat. Commun. 10, 3794 (2019).
  * Chen and Chen (2021) H.-B. Chen and Y.-N. Chen, “Canonical hamiltonian ensemble representation of dephasing dynamics and the impact of thermal fluctuations on quantum-to-classical transition,” Sci. Rep. 11, 10046 (2021).
  * Lin  _et al._ (2022) M.-C. Lin, P.-Y. Lo, F. Nori, and H.-B. Chen, “Precession-induced nonclassicality of the free induction decay of nv centers by a dynamical polarized nuclear spin bath,” J. Phys.: Condens. Matter 34, 505701 (2022).
  * Kropf  _et al._ (2016) C. M. Kropf, C. Gneiting, and A. Buchleitner, “Effective dynamics of disordered quantum systems,” Phys. Rev. X 6, 031023 (2016).
  * Chen  _et al._ (2018) H.-B. Chen, C. Gneiting, P.-Y. Lo, Y.-N. Chen, and F. Nori, “Simulating open quantum systems with hamiltonian ensembles and the nonclassicality of the dynamics,” Phys. Rev. Lett. 120, 030403 (2018).
  * Chen (2022) H.-B. Chen, “Effects of symmetry breaking of the structurally-disordered hamiltonian ensembles on the anisotropic decoherence of qubits,” Sci. Rep. 12, 2869 (2022).
  * Smithey  _et al._ (1993) D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, “Measurement of the wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum,” Phys. Rev. Lett. 70, 1244 (1993).
  * Breitenbach  _et al._ (1997) G. Breitenbach, S. Schiller, and J. Mlynek, “Measurement of the quantum states of squeezed light,” Nature 387, 471 (1997).
  * Lvovsky  _et al._ (2001) A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, “Quantum state reconstruction of the single-photon fock state,” Phys. Rev. Lett. 87, 050402 (2001).
  * Deléglise  _et al._ (2008) S. Deléglise, I. Dotsenko, C. Sayrin, J. Bernu, M. Brune, J.-M. Raimond, and S. Haroche, “Reconstruction of non-classical cavity field states with snapshots of their decoherence,” Nature 455, 510 (2008).
  * Flühmann and Home (2020) C. Flühmann and J. P. Home, “Direct characteristic-function tomography of quantum states of the trapped-ion motional oscillator,” Phys. Rev. Lett. 125, 043602 (2020).
  * Lutterbach and Davidovich (1997) L. G. Lutterbach and L. Davidovich, “Method for direct measurement of the wigner function in cavity qed and ion traps,” Phys. Rev. Lett. 78, 2547 (1997).
  * Banaszek  _et al._ (1999) K. Banaszek, C. Radzewicz, K. Wódkiewicz, and J. S. Krasiński, “Direct measurement of the wigner function by photon counting,” Phys. Rev. A 60, 674 (1999).
  * Bertet  _et al._ (2002) P. Bertet, A. Auffeves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, “Direct measurement of the wigner function of a one-photon fock state in a cavity,” Phys. Rev. Lett. 89, 200402 (2002).
  * Vlastakis  _et al._ (2013) B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, “Deterministically encoding quantum information using 100-photon schrödinger cat states,” Science 342, 607 (2013).
  * Winkelmann  _et al._ (2022) F.-R. Winkelmann, C. A. Weidner, G. Ramola, W. Alt, D. Meschede, and A. Alberti, “Direct measurement of the wigner function of atoms in an optical trap,” J. Phys. B: At. Mol. Opt. Phys. 55, 194004 (2022).
  * Ahmed  _et al._ (2021) S. Ahmed, C. Sánchez Muñoz, F. Nori, and A. F. Kockum, “Quantum state tomography with conditional generative adversarial networks,” Phys. Rev. Lett. 127, 140502 (2021).
  * Cuadras  _et al._ (2002) C. M. Cuadras, J. Fortiana, and J. A. Rodriguez-Lallena,  _Distributions with Given Marginals and Statistical Modelling_ (Springer, Netherlands, 2002).
  * Lin  _et al._ (2014) G. D. Lin, X. Dou, S. Kuriki, and J.-S. Huang, “Recent developments on the construction of bivariate distributions with fixed marginals,” J. Stat. Distrib. App. 1, 14 (2014).
  * Takeuchi (2010) T. T. Takeuchi, “Constructing a bivariate distribution function with given marginals and correlation: application to the galaxy luminosity function,” Mon. Not. Roy. Astron. Soc. 406, 1830 (2010).
  * Cicalese  _et al._ (2019) F. Cicalese, L. Gargano, and U. Vaccaro, “Minimum-entropy couplings and their applications,” IEEE Trans. Inf. Theory 65, 3436 (2019).
  * Ong  _et al._ (2021) S. H. Ong, R. C. Gupta, T. Ma, and S. Z. Sim, “Bivariate conway–maxwell poisson distributions with given marginals and correlation,” J. Stat. Theory Pract. 15, 10 (2021).
  * He  _et al._ (2016) K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in  _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (2016) pp. 770–778.
  * Xiong  _et al._ (2015) H.-N. Xiong, P.-Y. Lo, W.-M. Zhang, D. H. Feng, and F. Nori, “Non-markovian complexity in the quantum-to-classical transition,” Sci. Rep. 5, 13353 (2015).

## Acknowledgments

This work is supported by the National Science and Technology Council, Taiwan,
with Grants No. MOST 108-2112-M-006-020-MY2, MOST 109-2112-M-006-012, MOST
110-2112-M-006-012, MOST 111-2112-M-006-015-MY3, MOST 111-2112-M-A49-014, NSTC
112-2112-M-A49-019-MY3, NSTC 112-2123-M-006-001, MOST 112-2314-B-006-011, MOST
110-2224-E-007-003, and MOST 109-2222-E-006-005-MY2, partially by the Higher
Education Sprout Project, Ministry of Education to the Headquarters of
University Advancement at NCKU, and partially by the National Center for
Theoretical Sciences, Taiwan.

## Author contributions

H.-B.C. conceived and initiated the research. C.-H.L. developed the
methodology with input from B.-Y.T. under the supervision of C.-H.Y. and
H.-B.C. K.-L.L. optimized the model, collected the data, and analysed the
results under the supervision of C.-H.Y. and H.-B.C. P.-Y.L. and Y.-N.C
provided the physical meanings of the results. H.-B.C. drafted the manuscript
with input from P.-Y.L. and C.-H.Y. H.-B.C., Y.-N.C, and C.-H.Y. were
responsible for the integration among different research units.

## Competing interests

The authors declare no competing interests.