Iran University of Science and Technology] Department of Electrical
Engineering, Iran University of Science and Technology, Narmak, Tehran
16486-13114, Iran IR,NMR,UV

# A deep learning approach for inverse design of the metasurface for dual-
polarized waves

Fardin Ghorbani  Javad Shabanpour<EMAIL_ADDRESS>[  Sina
Beyraghi  Hossein Soleimani  Homayoon Oraizi  Mohammad Soleimani

###### Abstract

Compared to the conventional metasurface design, machine learning-based
methods have recently created an inspiring platform for an inverse realization
of the metasurfaces. Here, we have used the Deep Neural Network (DNN) for the
generation of desired output unit cell structures for both TE and TM polarized
waves which its working frequency can reach up to 45 GHz. To automatically
generate metasurfaces over wide frequencies, we deliberately design 8 annular
models; thus, each generated meta-atoms in our dataset can produce different
notches in our desired working frequency. Compared to the general approach,
whereby the final metasurface structure may be formed by any randomly
distributed “0” and “1”, we propose here a confined output configuration. By
confining the output, the number of calculations will be decreased and the
learning speed will be increased. Establishing a DNN-confined output
configuration based on the input data for both TE and TM polarized waves is
the novelty to generate the desired metasurface structure for dual orthogonal
polarizations. Moreover, we have demonstrated that our network can attain an
accuracy of 92%. Obtaining the final unit cell directly without any time-
consuming optimization algorithms for both TE and TM polarized waves, and high
average accuracy, open beneficial ways for the inverse metasurface design;
thus, the designer is required only to focus on the design goal.

###### keywords:

American Chemical Society, LaTeX

##  1 Introduction

Metamaterials have attracted widespread attention due to their peculiar assets
to modify the permittivity and permeability 1, 2, 27, 28, 29. Recently, many
novel functionalities have been implemented by metamaterials and their 2D
counterpart, metasurfaces, such as intelligent surfaces for communication3, 4,
real-time wavefront manipulation5, 6, 7, 8, perfect absorption9, 10, and
machine learning metasurface design11, 12.

However, all of these works are founded on conventional approaches including
trial-and-error methods, brute force optimization methods, and parameter
sweep, which are time-consuming processes. Therefore, to solve the above
challenges and to seek a fast, effective, and programmed route for designing a
metasurface, we have benefited from machine learning. Deep learning is an
efficient approach for learning the relationship between input and wanted
information from the samples of past experiences. To be more specific, deep
learning as a special section of machine learning can infer the basic rules
based on formerly specified data; then, for different assigned inputs, the
designed network can estimate reasonable decisions. With expanding development
of machine learning and its possible future applications to address some
important problems such as signal processing13 and through the wall imaging14,
we are now witnessing the opening of machine learning in wave-interaction
phenomena. Owing to its potential capacity to provide higher accuracy, less
design time, and enhance the productivity of a modeling procedure, machine
learning has been introduced in numerous electromagnetic phenomena, for
instance, all-dielectric metasurfaces15, antenna design16, 17, acoustic
metamaterials18, 19, and computational electromagnetics20, 21.

T. Cui et al. presented an inverse metasurface realization that can recognize
the internal principles between input metasurface structures and their
S-parameters with 76.5% of accuracy.22. Zhang et al. have proposed an approach
based on machine learning for designing anisotropic coding metasurfaces with a
connection between deep learning and BPSO to explore the ideal reflection
phases of two-unit cells for the desired target23. Shan et al. have introduced
conditional deep convolutional generative adversarial networks to code the
programmable metasurface for multiple beam steering, with accuracy higher than
94%.24. A machine learning-based inverse metasurface design has been provided,
in25 which is capable of directly computing the output metasurface structure
by entering the sought design targets into the model. Recently, a double deep
Q-learning network has been introduced for distinguishing the ideal type of
material and optimizing a hologram structure to increase its efficiency.26

Here, benefiting from Deep Neural Network (DNN), we propose a network scheme
for automatic metasurface design for dual-polarized waves with an average
accuracy of up to 92%. Our network is capable of generating any desired
metasurface based on the input data for both TE and TM polarized waves which
its working frequency can reach up to 45 GHz. The works presented in 22, 25
can generate the output structure in the frequency range of 5 to 20 GHz and
16-20 GHz, respectively. To broaden the working frequency, we consider 8
annular models with the purpose of generating single or multiple resonances in
the desired working frequency. Besides, to enhance the speed of the training,
shorten the number of computations, and boost the effectiveness of our
network, the output of the network is confined; thus, the DNN should generate
the metasurface structure by employing the proposed 8 annular models. We
demonstrate that the accuracy of this network also reaches 91%. Consuming less
computational resources, having a wide frequency band, generating desired
output metasurface without resorting to optimization procedure, and working
for both TE and TM polarized waves make our method promising for boosting the
speed of computations and designs.

##  2 Results

###  2.1 Metasurface Design

Figure 1 demonstrates the proposed 8 annular models, each of them is composed
of three layers of copper, a substrate (FR4 with permittivity of 4.2+0.025i,
and thickness of 1.5mm), and a ground plane to block the incident
electromagnetic waves. Each annular model is 1.6 mm and is composed of
$8\times 8$ lattices named as “1” and “0” to indicate the spaces with and
without the copper. Each meta-atom consists of $4\times 4$ arbitrarily
distributed annular models and the final unit-cell is 6.4 mm. Therefore, each
unit cell generated as an input of DNN comprises $32\times 32$ lattices with a
length of 0.2 mm. The reason for using 8 annular patterns is to produce single
or multiple resonances in the generated S-parameters in a broad frequency
range from 4 to 45 GHz. Since it is almost impossible to attain the
relationship between the input “0” and “1” matrix and its equivalent
S-parameter, the machine learning method can be considered as an encouraging
solution to shorten the computational operations for obtaining the ideal
results.

Figure 1: Diagram of the design procedure of confined- output configuration
and 8 annular models.

###  2.2 Deep learning

Nowadays, neural network algorithms have appeared to solve some fundamental
challenges especially in optimization and artificial intelligence. Schematic
representation of an artificial neuron is depicted in Figure 2, in which
${A_{n}}$ denotes the input neurons. Since each $A$ is connected to a weight,
the multiplication of $A_{i}$’s by $W_{i}$’s has emerged in the input of the
summation. By considering $\phi(x)$ as an activation function, eventually, the
output of this process is determined by:

$Y=\phi(\sum\limits_{i=1}^{n}W_{i}A_{i}+b_{i})$ (1)

In the above equation, $b$ indicates the bias value. Generally, artificial
neural networks are made of distinct layers of input, output, and a hidden
layer between them. By increasing the hidden layers, the complexity of the
network increases, and the neural network turns into a deep neural network.

Figure 2: A sketch of an artificial neuron.

###  2.3 Confined output configuration

In this paper, we have benefited from a DNN to find out the intrinsic
connections between the output generated metasurface and its S-parameters
features. To establish our dataset, by the means of the RAND function, a
collection of 2000 matrices that form the unit cells are created. Then, we use
CST Microwave Studio to determine the reflection characteristics of each unit
cell under the illumination of both TE and TM polarized waves. Then, by
linking CST MWS with MATLAB, the information regarding reflection
characteristics (resonance frequency, resonance depth, and resonance
bandwidth) are saved in a generated database. Note that for obtaining the
S-parameters of the unit cells in periodic structures, periodic boundary
conditions are adjusted along x- and y-directions, while open boundary
conditions are applied along the propagation of the incoming waves.

Table 1: Elaborate data of the confined output DNN configuration. Figure 3:
Two examples of reflection amplitude of confined output DNN configuration a,d)
final metasurface structure. b,e) simulated S-parameters under illumination of
TE and c,f) TM polarized waves. Figure 4: Two examples of reflection amplitude
of confined output DNN configuration a,d) final metasurface structure. b,e)
simulated S-parameters under illumination of TE and c,f) TM polarized waves.
Table 2: Pre-determined TE and TM input targets (TE / TM) for desired
S-parameters, that are demonstrated in Figure 3 and Figure 4.

Examples | resonance frequency (GHz) | resonance depth (dB) | resonance bandwidth (GHz)
---|---|---|---
Figure 3(a-c) | 6.5, 11.5, 24 / 6.5, 11.5 | -21.5, -22.5, -25 /-12, -21.5 | 0.3, 1, 0.6 / 0.2, 1
Figure 3(d-f) | 27.5 / 14 | -30 / -22 | 0.5 / 0.2
Figure 4(a-c) | 25.5, 40 / 38 | -12, -30 / -18 | 0.4, 0.6 / 0.5
Figure 4(d-f) | 24.5 / 10, 14.5, 33 | -34.5 / -18.5, -20, -14 | 0.5 / 0.2, 0.4, 0.3

Our dataset randomly generates 16 numbers from 1 to 8 to create $4\times 4$
matrix. Each number denotes one of the eight annular models; thus, the output
of our established model is a matrix of $32\times 32$. In the training step,
we have produced 2000 pairs of S-parameters (TE and TM) and their
corresponding metasurface structures (70% belongs to the training and the rest
of the data belongs to a testing set). Each generated unit cell can produce up
to eight resonances. Since we have extracted three features for each resonance
(resonance frequency, resonance depth, and resonance bandwidth), a vector with
a size of 24 has emerged in the input of our devised DNN. To enhance the
network speed, minimize the volume of computations, and boost the
effectiveness of our learning section, we have confined the output of the
network; so that, the final metasurface is formed employing our designed 8
annular models. Note that to create the final vector, eight annular models are
indicated by digital codes of “000” to “111”; thus, the output of the DNN
originates a vector with the length of 48. This confined output configuration
will reduce the volume of computations while maintaining the accuracy of the
network. The details of the confined output configuration are outlined in
Table I.

In our model, we have used dense and dropout layers successively as depicted
in Figure 1. To escape from the misdirecting in the learning section and
decreasing the chance of overfitting, we randomly neglected specified numbers
of neurons in the dropout layer. On the other hand, we have defined the dense
layer in such a way that each neuron has a connection to all the neurons in
the former layers. Throughout the training section, we have used Adam
optimizer to calculate the variances between the predetermined and output data
repeatedly. By defining a loss function as a Mean Square Error (Eq. 2), we
have observed the differences between the original and generated data; thus,
when it reaches the specified criterion, the training step will stop.

${\rm{MSE}}=\frac{1}{m}\sum\limits_{i=1}^{m}{{{({Y_{i}}-{{\hat{Y}}_{i}})}^{2}}}$
(2)

In the above equation, m, $Y_{i}$ and ${{{\hat{Y}}_{i}}}$ represent the number
of data points, observed and predicted values, respectively. To increase the
rate of accuracy, we adopted the sigmoid activation function in the layer
number of 11 (see Table 1) since our wanted output in the DNN is 0 or 1. Eq. 3
and Eq. 4 show the expression of the activation of relu and sigmoid operators.

$f(z)=\left\\{{\begin{array}[]{*{20}{l}}0&{for\,\,z<0}\\\ z&{for\,\,z\geq
0}\end{array}}\right.$ (3)

$\sigma(z)=\frac{1}{{1+{e^{-z}}}}$ (4)

When it comes to evaluating the model, different design targets of
S-parameters are expected to determine whether the designed network can
generate corresponding metasurface structures. To establish our model, the
Tensorflow and Keras frameworks have been adopted, while our network is
performed using Python 3.8.0. Finally, after the design steps are completed,
it is only needed to insert the wanted S-parameter features, and the designed
DNN can form the final unit cell in accordance with the learned information
obtained through the training process.

Table 3: Detailed information of training and evaluation time, in addition to
the model size for confined output DNN architecture.

| Confined output network
---|---
Training time | 19 minutes
Evaluation | 0.038 sec
Model size | 6 MB

To validate the efficiency of our presented DNN, four different examples are
provided in the following. The output matrix of the metasurfaces is generated
based on the input S-parameters data for both TE and TM polarized waves. The
designated reflection information, which is detailed in Table II, is as
follows for all the examples: [notches frequencies; notches depth; and notch
bandwidth]. Numbers before and after the slash symbol ( / ) demonstrate the
desired input data under the illumination of TE and TM-LP waves, respectively.
Then, by entering the final generated matrices into full-wave simulation
software, we obtained the simulated reflection amplitude under the
illumination of dual orthogonal polarizations. As an example, we intend to
design a metasurface with an S-parameter containing three resonances under TE
polarization while having two under -10 dB resonances under TM polarization
waves in specified frequencies (see the first row of Table II). Observe in
Figure 3(b,c) that the full-wave simulation successfully reaches the design
targets.

Figure 4: Diagrams of accuracy and mean square error according to 5000 Epochs
for confined output DNN configuration.

Similarly, for instance in the last example (see Figure 4(d-f)), a metasurface
is designed in such a way that its S-parameters contain one and three
resonances in pre-determined frequencies (see the last row of Table II). Full-
wave simulation results are in accordance with our sought design goals.
Moreover, the diagrams of the accuracy and loss function of our presented DNN
confined output architecture are depicted in Figure 5. In addition, the
details of training and evaluation time, and the model size for our proposed
DNN are presented in Table III. These results are captured utilizing Google
Colab with the model of Tesla T4 and 15 MB of RAM. As can be seen, the design
time of our approach to generating a unit cell is 0.038 sec. Therefore,
compared to the conventional method, which takes about 700 to 800 minutes, our
presented approach is much faster and more efficient. Note that we cannot
directly compare our method to other inverse designs of the metasurface since
the employed GPU and RAM are different.

Accordingly, we have proven that our machine learning-based approach is a
promising candidate for inverse design of the metasurfaces in terms of
calculation repetitions, accuracy and speed. Establishing a DNN-confined
output configuration based on the input data for both TE and TM polarized
waves is the innovation to form a specified metasurface structure for dual
orthogonal polarizations.

##  3 Conclusion

In conclusion, adopting a deep neural network, we have presented an inverse
metasurface design approach for dual orthogonal polarized waves. By merely
specifying four design goals for both TE and TM cases (number of notches,
notch frequencies, notch depths, and notch bandwidths), the designed DNN can
create the final metasurface structure in the output. To broaden the working
frequency, we have considered 8 annular models; so that the created unit cells
in the output can produce different resonances over wide frequencies up to 45
GHz. The numerical simulations illustrate that our DNN can successfully
generate the desired metasurface compared to our pre-determined designed
targets, with an average accuracy higher than 91%. We have shown that the
speed of our presented approach is much higher than conventional metasurface
design, and by proposing the confining output configuration, our approach
equips an encouraging platform as an efficient technique with respect to
computational repetitions, training and evaluation time, and average accuracy.
We believe that our used deep neural network approach is a suitable candidate
for inverse metasurface design for dual-polarized waves and complex wave-
interaction phenomena.

###  3.1 CONFLICT OF INTEREST

The authors declare that there is no conflict of interest

###  3.2 DATA AVAILABILITY

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

## References

  * 1 Rajabalipanah, H., Abdolali, A., Shabanpour, J., Momeni, A. & Cheldavi, A. Asymmetric spatial power dividers using phaseamplitude metasurfaces driven by huygens principle. ACS Omega 4, 14340–14352 (2019).
  * 2 Shabanpour, J. Full manipulation of the power intensity pattern in a large space-time digital metasurface: from arbitrary multibeam generation to harmonic beam steering scheme. Ann. Phys. 532, 2000321 (2020).
  * 3 Di Renzo, Marco, et al. ”Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead.” IEEE Journal on Selected Areas in Communications 38.11 (2020): 2450-2525.
  * 4 Di Renzo, Marco, et al. ”Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison.” IEEE Open Journal of the Communications Society 1 (2020): 798-807.
  * 5 Shabanpour, J. Programmable anisotropic digital metasurface for independent manipulation of dual-polarized THz waves based on a voltage-controlled phase transition of VO 2 microwires. J. Mater. Chem. 8, 7189–7199 (2020).
  * 6 Shabanpour, J., Beyraghi, S. & Cheldavi, A. Ultrafast reprogrammable multifunctional vanadium-dioxide-assisted metasurface for dynamic THz wavefront engineering. Sci. Rep. 10, 1–14 (2020).
  * 7 Javad Shabanpour, Sina Beyraghi, Fardin Ghorbani, and Homayoon Oraizi, ”Implementation of conformal digital metasurfaces for THz polarimetric sensing,” OSA Continuum 4, 1372-1380 (2021)
  * 8 Shabanpour, Javad, et al. ”Real-time multi-functional near-infrared wave manipulation with a 3-bit liquid crystal based coding metasurface.” Optics Express 29.10 (2021): 14525-14535.
  * 9 Landy, N. I., Sajuyigbe, S., Mock, J. J., Smith, D. R. & Padilla, W. J. Perfect metamaterial absorber. Phys. Rev. Lett. 100, 207402 (2008).
  * 10 Shabanpour, J., Beyraghi, S. & Oraizi, H. Reconfigurable honeycomb metamaterial absorber having incident angular stability. Sci. Rep. 10, 1–8 (2020).
  * 11 Gu, M., & Goi, E. Holography enabled by artificial intelligence. In Holography, Diffractive Optics, and Applications X (Vol. 11551, p. 1155102). International Society for Optics and Photonics (2020).
  * 12 Ghorbani, Fardin, et al. ”Deep neural network-based automatic metasurface design with a wide frequency range.” Scientific Reports 11.1 (2021): 1-8.
  * 13 Ghorbani, Fardin, et al. ”EEGsig machine learning-based toolbox for End-to-End EEG signal processing.” arXiv preprint arXiv:2010.12877 (2020).
  * 14 Ghorbani, Fardin, Hossein Soleimani, and Mohammad Soleimani. ”Deep Learning Approach for Target Locating in Through-the-Wall Radar under Electromagnetic Complex Wall.” arXiv preprint arXiv:2102.07990 (2021).
  * 15 An, S. et al. A deep learning approach for objective-driven all-dielectric metasurface design. ACS Photon. 6, 3196–3207 (2019).
  * 16 Cui, L., Zhang, Y., Zhang, R. & Liu, Q. H. A modified efficient KNN method for antenna optimization and design. IEEE Trans. Antennas Propag. 68, 6858–6866 (2020).
  * 17 Sharma, Y., Zhang, H. H. & Xin, H. Machine learning techniques for optimizing design of double T-shaped monopole antenna. IEEE Trans. Antennas Propag. 68, 5658–5663 (2020).
  * 18 Bacigalupo, Andrea, et al. ”Machine-learning techniques for the optimal design of acoustic metamaterials.” Journal of Optimization Theory and Applications (2019): 1-24.
  * 19 Wu, Rih-Teng, et al. ”Design of one-dimensional acoustic metamaterials using machine learning and cell concatenation.” Structural and Multidisciplinary Optimization (2021): 1-25.
  * 20 Yao, He Ming, et al. ”Machine learning methodology review for computational electromagnetics.” 2019 International Applied Computational Electromagnetics Society Symposium-China (ACES). Vol. 1. IEEE, 2019.
  * 21 Yao, He Ming, E. I. Wei, and Lijun Jiang. ”Two-step enhanced deep learning approach for electromagnetic inverse scattering problems.” IEEE Antennas and Wireless Propagation Letters 18.11 (2019): 2254-2258.
  * 22 Qiu, T. et al. Deep learning: a rapid and efficient route to automatic metasurface design. Adv. Sci. 6, 1900128 (2019).
  * 23 Zhang, Q. et al. Machine-learning designs of anisotropic digital coding metasurfaces. Adv. Theory Simul. 2, 1800132 (2019).
  * 24 Shan, T., Pan, X., Li, M., Xu, S. & Yang, F. Coding programmable metasurfaces based on deep learning techniques. IEEE J. Emerg. Sel. Topics Power Electron 10, 114–125 (2020).
  * 25 Shi, X., Qiu, T., Wang, J., Zhao, X. & Qu, S. Metasurface inverse design using machine learning approaches. J. Phys. D. 53, 275105 (2020).
  * 26 Sajedian, I., Lee, H. & Rho, J. Double-deep Q-learning to increase the efficiency of metasurface holograms. Sci. Rep. 9, 1–8 (2019).
  * 27 Rajabalipanah, Hamid, et al. ”Addition theorem revisiting for phase/amplitude-encoded metasurfaces: Asymmetric spatial power dividers.” arXiv preprint arXiv:1901.04063 (2019).
  * 28 Shabanpour, Javad, and Homayoon Oraizi. ”Some useful approximations for calculation of directivities of multibeam power patterns of large planar arrays.” arXiv preprint arXiv:2006.10423 (2020).
  * 29 Gholamian, Meysam, Javad Shabanpour, and Ahmad Cheldavi. ”Highly sensitive quarter-mode spoof localized plasmonic resonator for dual-detection RF microfluidic chemical sensor.” Journal of Physics D: Applied Physics 53.14 (2020): 145401.


# The role of tunneling in the ionization of atoms by ultrashort and intense
laser pulses

Gabriel M. Lando<EMAIL_ADDRESS>Université Paris-Saclay,
CNRS, LPTMS, 91405 Orsay, France

###### Abstract

Classically allowed transport is shown to compete with quantum tunneling
during the ionization of atoms by ultrashort and intense laser pulses, despite
Keldysh parameters smaller than unity. This is done by comparing exact
probability densities with the ones obtained from purely classical propagation
using the Truncated Wigner Approximation. Not only is classical transport
capable of moving trajectories away from the core, but it can also furnish
ionization probabilities of the same order as the quantum ones for intensities
currently employed in experiments. Our results have implications ranging from
a conceptual correction to semiclassical step models in strong-field physics
to the ongoing debate about tunneling time measurements in attoclock
experiments.

_Introduction –_ Tunneling is one of the most characteristically quantum
phenomena in nature. Often depicted as a near-magical violation of classical
conservation laws, it is the basic mechanism behind several physical and
chemical processes and technologies, such as: the stability of star cores [1],
the decay of radioactive elements [2], Josephson junctions [3], the transfer
of hydrogen atoms in chemical reactions [4, *hammes2006hydrogen, *Lowdin1963,
*Trixler2013], photosynthesis [8], electron transfer in proteins [9] and flash
memory cards [10, *Bez2003], and a multitude of others. It is then remarkable
that, to this day, the use of the term “tunneling” is not uniform among all
areas of physics and chemistry. This is especially true when the phenomenon
takes place in time domain, as opposed to the static tunneling tails
penetrating classically forbidden regions in quantum eigenstates [12].

For conservative systems with a single degree of freedom (DoF), it is tempting
to interpret time-dependent tunneling as the transport of portions of the wave
function through a tall potential barrier. However, this description does not
take momentum into account, and it is well-known that a more accurate picture
is given in terms of the Wigner function [13]. Here, it is possible to
directly visualize the full energy barrier in phase space, which forms a
separatrix, instead of its misleading projection in position space only. One
then sees that the Wigner function naturally extends over the momentum
direction and presents high-energy tails [13, 12, 14, 15]. If points are
sampled on these tails and evolved according to Hamilton’s equations, they
will emerge on the other side of the barrier through classically allowed
transport. This contrasts with the case of tunneling, where transport is
necessarily classically forbidden [16, 17, 18, 19].

The matter of whether classically allowed or forbidden transport is the
dominant mechanism for some particular phenomenon is relatively
straightforward for the simple systems in the above paragraph, but becomes a
lot more blurred in the high-dimensional or time-dependent case: The former is
deceptive because splitting phase space into disjoint regions is harder in
higher dimensions [20, 21]; The latter involves non-stationary potential
barriers through which energy is added to and/or subtracted from the system,
possibly triggering chaos even for a single DoF [22, 23, 24]. For both cases a
proper definition of tunneling is quite challenging, although progress has
been recently made in this direction [19].

The interactions between atoms and ultrashort, intense laser pulses belong to
the category of time-dependent processes with strongly chaotic classical
dynamics. Here, the mechanism referred to as tunneling ionization (TI) lies at
the heart of a plethora of phenomena, forming the first stage of several
semiclassical step models [25, 26]. There are good reasons to question the
role of tunneling in this context, as the overall classical effect of the
pulse is a chaotic shaking of classical trajectories, which are then perfectly
able to escape the core _via_ classically allowed ionization (CAI) [27]. The
number of trajectories scattered by CAI is even enough to reproduce higher-
harmonics generation (HHG) spectra semiclassically, despite the propagated
state being approximately the atom’s ground state [28]. Since semiclassical
calculations rely completely on classical trajectories, this would be
impossible if the transport pathways were truly dominated by TI [16].

In this manuscript, we do not make the standard distinction between vertical
(multi-photon) and horizontal (tunneling) ionization channels [29]. Instead,
we distinguish CAI from TI by direct comparisons with purely classical
simulations, where the distribution of initial trajectories does not come from
any tweaking [30, 31, 32], but from the system’s _true_ ground state. Since
tunneling is stricter and easier to spot in systems with a single DoF, our
simulations are performed using the improved soft-core potential of Majorosi
_et al_ [33] coupled to an intense, ultrashort and linearly polarized (LP)
laser pulse in the near-infrared range. Contrary to intuition, but in line
with [20, 28], we demonstrate that CAI and TI are deeply intertwined, and
possibly inseparable. Among the implications of these results lie a correction
to semiclassical step models, since ionization is mostly unrelated to
tunneling, and an added difficulty to the ongoing debate on “tunneling” time
measurements [34, 35, 30, 36, 37, 32, 38, 39, 40].

_Model system–_ One-dimensional systems allow for visualization ease, but
often disagree with full 3 DoF simulations – at least quantitatively. The
popular soft-core potential, for example, has been shown to largely
overestimate ground-state depletion and other expectation values [33].
Nevertheless, recent work by Majorosi _et al_ has shown that the improved
soft-core system (in atomic units)

$H_{\rm{iSC}}(p_{x},x)=p_{x}^{2}/2-Z\left(Z^{-2}+4x^{2}\right)^{-1/2}\quad$
(1)

offers strikingly accurate estimates when compared to simulations in full
dimensionality using the exact Coulomb potential [33]. We therefore adopt this
potential for hydrogen, setting $Z=1$, and couple it to a linearly polarized
laser pulse in the dipole approximation:

$V(x;t)=x\,\sqrt{I_{0}}f(t)\sin\omega t\quad,$ (2)

where the envelope is given by $f(t)=\exp-((t-t_{0})/\tau)^{2}$. The optical
frequency of the pulse is taken as $\omega\approx
0.057\,\text{au}=780\,\text{nm}$, lying in the near-infrared. The FWHM is
$\tau\approx 4\,\text{fs}$, such that centering the pulse at $t_{0}\approx
10\,\text{fs}$ covers it completely for $t\in[0,20]\,\text{fs}$. Intensities
$I_{0}$ will vary between $1$ and $4\times 10^{14}\,\text{W/cm}^{2}$. The
initial state is chosen as the exact ground state for (1), which we obtain
numerically by diagonalizing the hamiltonian on a position grid ranging from
$-1000$ to $1000\,\text{au}$, with step-size $\Delta x\approx 0.1\,\text{au}$.

We note that the substitution of (1) by any other typical atomic potential, as
well as performing simulations with more DoF, do not change the message
contained in this manuscript: We are concerned only with the differences
between classical and quantum, and dynamics due to LP pulses is essentially
restricted to a single DoF [33]. Another important aspect is that quantum-
classical agreement for systems that undergo tunneling is harder to achieve in
one DoF than in more DoF. This was already remarked for a static electric
field in [20], where the author correctly states that trajectory “leakage”,
which provides the trajectories necessary to reproduce ionization
semiclassically, is more abundant in two DoF than in one. Thus, our choice of
(1) as a model both simplifies interpretation and _overestimates_ the role
played by TI.

_Methods–_ Quantum time-evolution is performed by solving the time-dependent
Schrödinger equation, where a split-operator method [41] using Blanes and
Moan’s 6th order algorithm is employed [42]. The time-step chosen is $\Delta
t=0.2\,\text{au}$, corresponding to $1/550$ of an optical cycle. To mitigate
wave reflections, we use a $\cos^{1/8}$ boundary mask placed at 10% of the
grid’s extremities [43, 44].

Quantum dynamics in phase space, _i.e._ the time evolution of the Wigner
function

$W_{0}(p_{x},x)=\frac{1}{\pi\hbar}\int_{\mathbb{R}}\text{d}\gamma\,\langle
x-\gamma|\psi_{0}\rangle\langle\psi_{0}|x+\gamma\rangle e^{2i\gamma
p/\hbar}\quad,$ (3)

where $|\psi_{0}\rangle$ is the system’s ground state, is dictated by Moyal’s
equation [45, 46, 15, 24]. Here, unlike in other formulations, the $\hbar\to
0$ limit of time-evolution is well-defined and results in a von Neumann-like
equation, with solution given by the Truncated Wigner Approximation (TWA)

$w(p_{x},x;t)=(W_{0}\circ\varrho_{-t})(p_{x},x)\quad.$ (4)

In the above, $\varrho_{-t}$ is the hamiltonian flow bringing the point
$(p_{x},x)$ at $\tau=0$ to its value at $\tau=-t$. Computing this backwards
flow requires us to employ a time-mirrored laser pulse in Hamilton’s
equations, which we solve using an adaptive-step algorithm with automatic
stiffness detection [47, 44].

Expectation values using the TWA follow the recipe of the Wigner formalism:

$\langle
A(t)\rangle_{\rm{classical}}=\int_{\mathbb{R}^{2}}\text{d}x\,\text{d}p_{x}\,w(p_{x},x;t)A(p_{x},x)\quad,$
(5)

where $A$ is the Weyl transform of operator $\hat{A}$ [15]. We note that the
usual way of computing classical expectation values is to transfer the time-
dependence from $w$ to $A$, and then propagate forward in time [48]. We go
through the trouble of negative times because, just as the marginals of Wigner
functions describe quantum probability densities (PDs), the marginals obtained
from the TWA can be seen as their classical equivalents. Comparisons between
quantum and classical PDs will be our tool to determine whether or not TI is
dominant: If this is the case, the classical PDs will have empty regions when
compared to the quantum ones, indicating that the corresponding phase-space
domain cannot be accessed by classical trajectories; If not, quantum and
classical PDs will be non-zero on the same domain.

In practice, since (5) is performed by Monte Carlo, the classical (position)
PDs are nothing but histograms of final positions in the TWA. These, in turn,
come from a classically evolved, non-interacting gas of initial points usually
chosen by importance sampling. Since we are interested in evolving the ground
state Wigner function $W_{0}$, sampling the initial points according to
$|W_{0}|$ is very efficient, and we do so by employing a simple Metropolis-
Hastings scheme [49, 50]. The absolute value is necessary because the ground
state of (1) is not a gaussian and, therefore, $W_{0}$ presents negative-
valued regions (see Fig. 1 ahead) [51]. It should perhaps be mentioned that
this association between a quantum state and a classical ensemble of points is
a sensible way to establish quantum-classical correspondence for dynamics,
being far more general than Ehrenfest’s theorem [52, 53, 54, 19].

Figure 1: The ground state Wigner function, $W_{0}$. The solid black line is
the zero energy contour, so orbits lying inside (outside of) it are bounded
(scattered). The arrows mark orbits of type N1 (magenta), N2 (green) and P
(orange).

_Simulations–_ In the case of potentials such as (1), which is composed of
bounded and scattering regions (corresponding to discrete and continuous
quantum energy spectra), the initial points can have both closed and open
trajectories. For (1), closed and open trajectories have negative and positive
energies, respectively, and we shall refer to them as being of type N or type
P.

Type N trajectories are all periodic. In the absence of the pulse some of
them, which we call type N1, naturally extend far beyond the mean ground state
radius, orbiting the origin with large periods. Other type N trajectories,
which we shall call N2, remain near the origin and have small orbital periods.
Since the laser pulse is in the near-infrared range, it acts on type N2 and N1
trajectories adiabatically and non-adiabatically, respectively. Thus,
elementary classical perturbation theory [55, 23] tells us that the pulse will
keep type N2 trajectories mostly untouched, and scatter away some of type N1.
The former mechanism is responsible for the classical preservation of the
ground state, since it conserves trajectories near the origin, and the latter
is the one behind CAI. See Fig. 1 for a visual depiction of these concepts.

Figure 2: (a) Final Wigner function and (b) TWA after a pulse with intensity
$I_{0}=2.0\times 10^{14}\,\text{W/cm}^{2}$. Blue (red) colors in the Wigner
function are positive (negative), and the TWA is built on around $10^{6}$
trajectories. The trajectories in panel (b) were all ionized by the laser
field, since the ones lying on the tails of the Wigner function are absent in
the set of initial points (see text). Note that the Wigner function and the
TWA have essentially the same spread.

The situation is very different with type P trajectories, whose initial points
lie on the Wigner function’s “tails” (see Fig. 1). Due to the exponential
decay of the tails, less than 0.01% of the initial points form trajectories of
type P. As these trajectories are not bound to the core, they are scattered
away in a matter of a few attoseconds, even in the absence of the laser pulse.
Moreover, by galilean invariance, we could simply translate the pulse to a
time where type P trajectories would have already diverged, so there are many
reasons to consider that they are irrelevant. After verifying that our results
are unchanged whether or not we include them, we simply remove them. This
filtering allows us to track classical ionization probabilities exclusively to
the trajectories that _truly_ ionized, _i.e._ they started bounded and were
scattered by the field. In Fig. 2 we show the Wigner function at the end of a
pulse with intensity $I_{0}=2.0\times 10^{14}\,\text{W/cm}^{2}$ together with
the corresponding TWA, calculated using a filtered set of around $10^{6}$
initial points.

Figure 3: (Left) Heatmap of the probability of finding the electron as a
function of space and time during the laser pulse. The color scale in the left
panel is logarithmic, with green equal to $10^{-2}$ and purple to $10^{-5}$.
(Right) Classical (orange) and quantum (black) position PDs at the end of the
laser pulse. Note that the peak probability is close to $0.5$, but ionization
probabilities are no larger than $0.0006$, requiring a massive zoom in order
to be seen.

Fig. 2 shows that CAI is taking place, since the TWA extends far beyond the
core, and we now move on to quantify how much this classical spread accounts
for total ionization. In the left panel of Fig. 3 we display the quantum PD in
the system as a function of time for the same intensity as Fig. 2. Results are
shown from $5\,\text{fs}$ (which is when the laser pulse starts visibly acting
on the state) to $20\,\text{fs}$. In the right panel we show the final quantum
PD, $P(x)$, together with its classical approximation, obtained from the TWA’s
position marginal.

_Discussion–_ In most situations where purely classical approximations are
employed, what is expected is only a very rough sketch of what is obtained
quantum mechanically. Classical physics cannot reproduce quantum
superposition, which plays a fundamental role in several processes in strong-
field physics (HHG, for instance, relies strongly on it [31]). This makes it
impossible to reconstruct parts of PD that are strongly dependent on phase
interference. The results of Fig. 3 are surprising because the classical
approximation is far more than a rough sketch: Not only does it correctly
estimate the (minuscule) order of atomic ionization probability, but it also
highlights the overall structure of the quantum result. Most importantly, the
classical PD extends over essentially the same range as the quantum one
despite a Keldysh parameter smaller than unity, namely
$\gamma_{\text{K}}\approx 0.78$, which should be indicative of a strong
presence of tunneling [29].

The sharp peaks displayed by the classical PD in the right panel of Fig. 3 are
due to the “whorls” and “tendrils” of Fig. 2, _i.e._ the initial state is
deformed into filaments that are sheared and folded, and the classical
position marginals have peaks on the filaments lying perpendicularly to the
momentum axis [56, 57]. This filamentary structure is usually not enough to
reproduce peak heights, and one must resort to semiclassical approximations:
Quantum interference is then achieved by a rigorous endowing of classical
trajectories with accumulated phases, which superpose and correct peak
intensities [58, 59]. What is fundamental to our purposes is that only the
extremities of the PD in Fig. 3, which are barely visible, are not classically
accessible – These are the ones that emerge exclusively from TI.

Figure 4: Quantum (black) and classical (orange) PDs at the end of the laser
pulse for: (a) $I_{0}=1.0\times 10^{14}\,\text{W/cm}^{2}$; (b)
$I_{0}=4.0\times 10^{14}\,\text{W/cm}^{2}$. Classical calculations employed
$10^{6}$ trajectories. The arrows in panel (b) point at regions where TI can
be safely disentangled from CAI.

The success of both classical and semiclassical mechanics, however, depends on
the existence of available classical trajectories. After they are all
scattered away, TI becomes the only available escape pathway from the core, as
can be seen for static fields or long laser pulses [20]. If the pulses are
strong, trajectories will be scattered faster, such that very strong fields
should also pose problems for classical/semiclassical propagation 111Note that
this is in line with standard Keldysh theory, but for different reasons..
Thus, in Fig. 4 we display the final probability densities for two different
laser pulse intensities: One weaker ($\gamma_{\text{K}}\approx 1.1$) and the
other one stronger ($\gamma_{\text{K}}\approx 0.54$) than in Figs. 2 and 3.

As we can see in Fig. 4(a), for “weak” fields the quantum PD is completely
supported by the classical one, which is almost one order of magnitude larger.
This shows that, in the limit of weak fields, CAI is notoriously dominant over
TI, so much so that it even resembles a case of strong localization [61].
Whether or not semiclassical corrections will be able to suppress the
classical PD and reproduce the peaks depends on how far one is in the
semiclassical regime, _i.e._ how large the actions associated to the
trajectories are with respect to $\hbar$, but a perfect reproduction of
probability densities is not the objective of this manuscript. What is
important is that there’s no region in position space that classical
trajectories didn’t reach, so TI cannot be dominant. In Fig. 4(b) however, we
note the increasing presence of peaks falling outside the reach of classical
trajectories, and TI is seen to arise from outside in. Using Fig. 3, we can
describe this quite easily: Some portions of the wave function that ionize
near the maxima of the laser pulse, with the absolute maximum at
$t_{0}=10\,\text{fs}$, go farther. They provide peaks at the PD’s extremities
that cannot be fully reached by CAI, although this is hard to see for the
intensity used in Fig. 3 and much clearer in Fig. 4(b). In addition to the
classically unreachable extremities, it is also unlikely that semiclassical
approximations will be able to resolve most of the peaks in Fig. 4(b), as it
is clear that both TI and CAI are taking place.

Interestingly, if we interpret Fig. 4(b) together with Fig. 3, we see that the
purest tunneling contributions can be tracked to portions that ionize and
never recombine – In the terminology of HHG literature, the classical
trajectories near these portions (if there are any) never _recollide_ with the
core [62, 29, 63]. The recolliding trajectories abide, as expected, within a
ball of quiver radius $\alpha_{0}=\sqrt{I_{0}}/\omega^{2}$ [29]. As can be
seen in the Figs. 2, 3 and 4, however, this region is always well-populated by
classical trajectories. This explains why HHG spectra could be obtained in
[28] through semiclassical approximations: Tunneling contributions might even
be present near the core, but the classical trajectories in the region were
enough to reproduce spectra exclusively from CAI.

It is then reasonable that the first step in semiclassical step models be
called “ionization” instead of “tunneling ionization”, since one cannot
effectively know whether it came from tunneling or not. This, moreover, adds
to the previously raised concerns about the difficulty of theoretically and
experimentally resolving tunneling time ambiguities, since the contributions
from TI and CAI might be impossible to disentangle. In the end, it is possible
that the measured angular delays in attoclocks, often interpreted as the “time
spent inside a barrier”, are predominantly due to CAI and directly traceable
to transport properties of classical trajectories. This would render
“tunneling times” essentially unrelated to tunneling.

_Conclusions–_ We have used a one-dimensional atom to bring forward the fact
that isolating tunneling from classically allowed ionization can be extremely
hard, even for the simplest of systems, when they are interacting with intense
and short laser pulses. Our results show that what is known as “tunneling
ionization” in semiclassical step models in atomic physics has a major
classical fingerprint, with direct consequences to the proper interpretation
of higher-harmonics generation and tunneling times, among others.

_Acknowledgements–_ I thank Alfredo Ozorio de Almeida, Andrew Hunter, Denis
Ullmo, Frank Großmann, Jessica Almeida, Jonathan Dubois, Olivier Giraud, Peter
Schlagheck, Sebastian Gemsheim and Steven Tomsovic for many stimulating
discussions. I also thank Jan-Michael Rost and the hospitality of the Max
Planck Institute for the Physics of Complex Systems, where the initial stages
of this work were carried out.

## References

  * Itoh  _et al._ [1979] N. Itoh, H. Totsuji, S. Ichimaru, and H. E. Dewitt, Enhancement of thermonuclear reaction rate due to strong screening. II - Ionic mixtures, Astrophys. J. 234, 1079 (1979).
  * Gurney and Condon [1928] R. W. Gurney and E. U. Condon, Wave mechanics and radioactive disintegration, Nature 122, 439 (1928).
  * Josephson [1962] B. Josephson, Possible new effects in superconductive tunnelling, Physics Letters 1, 251 (1962).
  * Kohen and Klinman [1999] A. Kohen and J. P. Klinman, Hydrogen tunneling in biology, Chemistry & biology 6, R191 (1999).
  * Hammes-Schiffer [2006] S. Hammes-Schiffer, Hydrogen tunneling and protein motion in enzyme reactions, Accounts of chemical research 39, 93 (2006).
  * Löwdin [1963] P.-O. Löwdin, Proton tunneling in dna and its biological implications, Reviews of Modern Physics 35, 724 (1963).
  * Trixler [2013] F. Trixler, Quantum tunnelling to the origin and evolution of life, Curr. Org. Chem. 17, 1758 (2013).
  * Peters  _et al._ [1978] K. Peters, P. Avouris, and P. Rentzepis, Picosecond dynamics of primary electron-transfer processes in bacterial photosynthesis, Biophysical Journal 23, 207 (1978).
  * Gray and Winkler [2003] H. B. Gray and J. R. Winkler, Electron tunneling through proteins, Quarterly reviews of biophysics 36, 341 (2003).
  * Esaki [1974] L. Esaki, Long journey into tunneling, Science 183, 1149 (1974).
  * Bez  _et al._ [2003] R. Bez, E. Camerlenghi, A. Modelli, and A. Visconti, Introduction to flash memory, Proceedings of the IEEE 91, 489 (2003).
  * Berry and Mount [1972] M. V. Berry and K. Mount, Semiclassical approximations in wave mechanics, Reports on Progress in Physics 35, 315 (1972).
  * Balazs and Voros [1990] N. Balazs and A. Voros, Wigner’s function and tunneling, Annals of Physics 199, 123 (1990).
  * Berry [1977] M. V. Berry, Semi-classical mechanics in phase space: a study of wigner’s function, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 287, 237 (1977).
  * De Almeida [1998] A. M. O. De Almeida, The weyl representation in classical and quantum mechanics, Physics reports 295, 265 (1998).
  * Maitra and Heller [1997] N. Maitra and E. Heller, Barrier tunneling and reflection in the time and energy domains: The battle of the exponentials, Physical review letters 78, 3035 (1997).
  * Grossmann and Heller [1995] F. Grossmann and E. J. Heller, A semiclassical correlation function approach to barrier tunneling, Chemical physics letters 241, 45 (1995).
  * Grossmann [2000] F. Grossmann, Semiclassical real-time tunneling by multiple spawning of classical trajectories, Physical review letters 85, 903 (2000).
  * Wang and Tomsovic [2021] H. Wang and S. Tomsovic, Semiclassical propagation of coherent states and wave packets: hidden saddles, arXiv preprint arXiv:2107.08799 (2021).
  * Spanner [2003] M. Spanner, Strong field tunnel ionization by real-valued classical trajectories, Physical review letters 90, 233005 (2003).
  * Zagoya  _et al._ [2014a] C. Zagoya, L. S. Schulman, and F. Grossmann, Interference nature of quantum breather oscillation, Journal of Physics A: Mathematical and Theoretical 47, 165102 (2014a).
  * Chirikov [1979a] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Physics reports 52, 263 (1979a).
  * Arnold [1989] V. I. Arnold,  _Mathematical Methods of Classical Mechanics_ , 2nd ed. (Springer, 1989).
  * de Almeida [1990] A. M. O. de Almeida,  _Hamiltonian Systems: Chaos and Quantization_ , 2nd ed. (Cambridge University Press, 1990).
  * Corkum [1993] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Physical review letters 71, 1994 (1993).
  * Lewenstein  _et al._ [1994] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’huillier, and P. B. Corkum, Theory of high-harmonic generation by low-frequency laser fields, Physical Review A 49, 2117 (1994).
  * Dubois [2019] J. Dubois,  _Electron dynamics for atoms driven by intense and elliptically polarized laser pulses_ , Ph.D. thesis, Aix-Marseille (2019).
  * Zagoya  _et al._ [2014b] C. Zagoya, J. Wu, M. Ronto, D. Shalashilin, and C. F. de Morisson Faria, Quantum and semiclassical phase-space dynamics of a wave packet in strong fields using initial-value representations, New Journal of Physics 16, 103040 (2014b).
  * Großmann [2013] F. Großmann,  _Theoretical Femtosecond Physics: Atoms and Molecules in Strong Laser Fields_ , 2nd ed. (Springer, 2013).
  * Ni  _et al._ [2016] H. Ni, U. Saalmann, and J.-M. Rost, Tunneling ionization time resolved by backpropagation, Physical review letters 117, 023002 (2016).
  * van de Sand and Rost [1999] G. van de Sand and J. M. Rost, Irregular orbits generate higher harmonics, Physical review letters 83, 524 (1999).
  * Hofmann  _et al._ [2019] C. Hofmann, A. S. Landsman, and U. Keller, Attoclock revisited on electron tunnelling time, Journal of Modern Optics 66, 1052 (2019).
  * Majorosi  _et al._ [2018] S. Majorosi, M. G. Benedict, and A. Czirják, Improved one-dimensional model potentials for strong-field simulations, Physical Review A 98, 023401 (2018).
  * Landsman  _et al._ [2014] A. S. Landsman, M. Weger, J. Maurer, R. Boge, A. Ludwig, S. Heuser, C. Cirelli, L. Gallmann, and U. Keller, Ultrafast resolution of tunneling delay time, Optica 1, 343 (2014).
  * Torlina  _et al._ [2015] L. Torlina, F. Morales, J. Kaushal, I. Ivanov, A. Kheifets, A. Zielinski, A. Scrinzi, H. G. Muller, S. Sukiasyan, M. Ivanov, _et al._ , Interpreting attoclock measurements of tunnelling times, Nature Physics 11, 503 (2015).
  * Pollak [2017] E. Pollak, Quantum tunneling: the longer the path, the less time it takes, The journal of physical chemistry letters 8, 352 (2017).
  * Rost and Saalmann [2019] J. M. Rost and U. Saalmann, Attoclock and tunnelling time, Nature Photonics 13, 439 (2019).
  * Kheifets [2020] A. S. Kheifets, The attoclock and the tunneling time debate, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 072001 (2020).
  * Sainadh  _et al._ [2020] U. S. Sainadh, R. Sang, and I. Litvinyuk, Attoclock and the quest for tunnelling time in strong-field physics, Journal of Physics: Photonics 2, 042002 (2020).
  * Hofmann  _et al._ [2021] C. Hofmann, A. Bray, W. Koch, H. Ni, and N. I. Shvetsov-Shilovski, Quantum battles in attoscience: tunnelling, The European Physical Journal D 75, 1 (2021).
  * Feit  _et al._ [1982] M. Feit, J. Fleck Jr, and A. Steiger, Solution of the schrödinger equation by a spectral method, Journal of Computational Physics 47, 412 (1982).
  * Blanes and Moan [2002] S. Blanes and P. C. Moan, Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods, Journal of Computational and Applied Mathematics 142, 313 (2002).
  * Lorin  _et al._ [2009] E. Lorin, S. Chelkowski, and A. Bandrauk, Mathematical modeling of boundary conditions for laser-molecule time-dependent schrödinger equations and some aspects of their numerical computation—one-dimensional case, Numerical Methods for Partial Differential Equations: An International Journal 25, 110 (2009).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM review 59, 65 (2017).
  * Groenewold [1946] H. J. Groenewold, On the principles of elementary quantum mechanics, in  _On the principles of elementary quantum mechanics_ (Springer, 1946) pp. 1–56.
  * Moyal [1949] J. E. Moyal, Quantum mechanics as a statistical theory, in  _Mathematical Proceedings of the Cambridge Philosophical Society_ , Vol. 45 (Cambridge University Press, 1949) pp. 99–124.
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, Differentialequations.jl – a performant and feature-rich ecosystem for solving differential equations in julia, The Journal of Open Research Software 5 (2017), exported from https://app.dimensions.ai on 2019/05/05.
  * Mittal  _et al._ [2020] K. M. Mittal, O. Giraud, and D. Ullmo, Semiclassical evaluation of expectation values, Physical Review E 102, 042211 (2020).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953).
  * Hastings [1970] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 236, 97 (1970).
  * Hudson [1974] R. L. Hudson, When is the wigner quasi-probability density non-negative?, Reports on Mathematical Physics 6, 249 (1974).
  * Ballentine  _et al._ [1994] L. E. Ballentine, Y. Yang, and J. Zibin, Inadequacy of ehrenfest’s theorem to characterize the classical regime, Physical review A 50, 2854 (1994).
  * Drobnỳ  _et al._ [1997] G. Drobnỳ, A. Bandilla, and I. Jex, Quantum description of nonlinearly interacting oscillators via classical trajectories, Physical Review A 55, 78 (1997).
  * Lasser and Lubich [2020] C. Lasser and C. Lubich, Computing quantum dynamics in the semiclassical regime, Acta Numerica 29, 229 (2020).
  * Chirikov [1979b] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Phys. Rep. 52, 263 (1979b).
  * Berry [1979] M. V. Berry, Evolution of semiclassical quantum states in phase space, Journal of Physics A: Mathematical and General 12, 625 (1979).
  * Berry  _et al._ [1979] M. V. Berry, N. L. Balazs, M. Tabor, and A. Voros, Quantum maps, Annals of Physics 122, 26 (1979).
  * Maslov and Fedoriuk [1981] V. P. Maslov and M. V. Fedoriuk,  _Semi-Classical Approximation in Quantum Mechanics_ (Springer, 1981).
  * Lando  _et al._ [2019] G. M. Lando, R. O. Vallejos, G.-L. Ingold, and A. M. O. de Almeida, Quantum revival patterns from classical phase-space trajectories, Physical Review A 99, 042125 (2019).
  * Note [1] Note that this is in line with standard Keldysh theory, but for different reasons.
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Physical review 109, 1492 (1958).
  * Protopapas  _et al._ [1996] M. Protopapas, D. Lappas, C. H. Keitel, and P. L. Knight, Recollisions, bremsstrahlung, and attosecond pulses from intense laser fields, Physical Review A 53, R2933 (1996).
  * Dubois  _et al._ [2020] J. Dubois, C. Chandre, and T. Uzer, Envelope-driven recollisions triggered by an elliptically polarized pulse, Physical Review Letters 124, 253203 (2020).


\intertext{As before, the only candidate for a joint distribution with finite score is $\delta_x(X) e(Z \mid X)$. Note that the marginal on $Z$ for this distribution is itself, since $\int_x \delta_x(X) e(Z \mid X)\;\mathrm dx = e(Z \mid x)$. Thus, our equation becomes}
&= \Ex_{\delta_x(X) e(Z \mid X)} \left[ \beta \log \frac{e(Z \mid x)}{p(z)} + \log \frac{\delta_x(X) e(Z \mid X)}{e(Z \mid x) d(x \mid Z)} \right] \\
&= \Ex_{e(Z \mid x)} \left[ \beta \log \frac{e(Z \mid x)}{p(Z)} + \log \frac{1}{ d(x \mid Z)} \right]\\
&= \kldiv{e(Z|\,x)}{p} + \mathrm{Rec}_{e,d}(x) \\
&= -\beta\text{-}\mathrm{ELBO}_{p,e,d}(x).
\end{align*}

In the main text, we defined $\PDGof{\Psi}$ to be the PDG with edges $\{ \raisebox{-0.3ex}{$\smash{\stackrel{J}{\rightarrow}}$} \mathbf X_J \}_{\mathcal J}$, cpds $p_J(\mathbf X_J) \propto \phi_J(\mathbf X_J)$, and weights $\alpha_J, \beta_J := \theta_J$.
Let $\theelt(\{x\}) := x$ be a function that extracts the unique element singleton set.

It was shown by richardson2020probabilistic (Corolary 4.4.1) that

\[ \theelt \bbr{\dg M_\Psi}^*_1 = \Pr\nolimits_{\Phi, \theta}(\mat x)
= \frac{1}{Z_\Psi} \prod_{J} \phi_J(\mat x_J)^{\theta_J}. \]

Recall the statement of Prop 4.6 from richardson2020probabilistic:

\begin{equation}\label{eqn:nice-score-repeated}
\bbr{\dg M}_\gamma(\mu) = \Ex_{\mat w \sim \mu}\! \Bigg\{ \sum_{ X \xrightarrow{\!\!L} Y  } \bigg[\,
\!\beta_L \log \frac{1}{\bp(y^{\mat w} |x^{\mat w})} +
{\color{blue!50!red!90!black}(\gamma\alpha_L - \beta_L ) \log \frac{1}{\mu(y^{\mat w} |x^{\mat w})}} \bigg] -
\gamma \log \frac{1}{\mu(\mat w)}  \Bigg\}, \\
\end{equation}

where $x^{\mat w}$
and $y^{\mat w}$ are the respective values of the variables $X$ and $Y$ in the world $\mat w$.
Note that if  $\gamma = 1$, and $\alpha,\beta$ are both equal to $\theta$ in $\PDGof{\Psi}$,

the  middle term (in purple) is zero. So in our case, since the edges are $\{ \xrightarrow{J} \mathbf X_J \}$ and $\bp[J](\mat X_J) = \phi_J(\mathbf X_J)$, (<ref>) reduces to the standard variational free energy

\begin{align*}
\VFE_\Psi(\mu)
&= \Ex_{\mu} \left[ ~\sum_{J\in \mathcal J} \theta_J \log  \frac1{\phi_J(\mat X_J)}\right]-\H(\mu) \numberthis\label{eq:vfe}\\
\qquad
\Ex_{\mu}
% \sum_{J\in \mathcal J} {\varphi_J(\mat X_J)}
\big\langle \boldsymbol\varphi,\, \boldsymbol\theta \big\rangle_{\mathcal J}
- \H(\mu),
\quad\text{where}~\varphi_J(\mat X_J) := \log \frac1{\phi_J(\mat X_J)}.
\end{align*}

By construction, $\Pr_\Psi$ uniquely minimizes $\VFE$.
The 1-inconsistency, $\aar{\dg M_\Psi}$ is the minimum value attained. We calculate:

\begin{align*}
% \aar{(\UPDGof{\Phi}, \theta, \theta)}_1
\aar{\dg M}_1
% &=\bbr{(\UPDGof{\Phi}, \theta, \theta)}_1\Big(\Pr\nolimits_{\Phi, \theta}(\mat w) \Big)\\
&= \VFE_\Psi(\Pr\nolimits_\Psi) \\
%\frac{1}{Z_\Phi} \prod_j \phi_j(\mat w_j)^{\theta_j}
% \Ex_{\mat w \sim \mu}\! \Bigg\{ \sum_{X \xrightarrow{\!\!L} Y} \bigg[\,
% \!\beta_L \log \frac{1}{\bp(y^{\mat w} |x^{\mat w})}
% \bigg] - \log \frac{1}{\Pr\nolimits_{\Phi, \theta}(\mat w) }  \Bigg\}
% & \Big[	 ~\text{by \eqref{eqn:nice-score-repeated}}~	\Big]\\
% \Ex_{\mat w \sim \mu}\! \Bigg\{ \sum_{X \xrightarrow{\!\!L} Y} \bigg[\,
\Ex_{\mat x \sim \mu}\! \Bigg\{ \sum_{J \in \mathcal J} \bigg[\,
\!\theta_J \log \frac{1}{\phi_J(\mat x_J)}
% (\alpha_L - \beta_L ) \log \frac{1}{\mu(y^{\mat w} |x^{\mat w})}
\bigg] - \log \frac{1}{\Pr\nolimits_{\Phi, \theta}(\mat x) }  \Bigg\}
& \Big[	 ~\text{by \eqref{eq:vfe}}~	\Big]\\
\Ex_{\mat x \sim \mu}\! \Bigg\{ \sum_{J\in \mathcal J} \bigg[\,
\!\theta_J \log \frac{1}{\phi_J(\mat x_J)}
\bigg] - \log \frac{Z_\Psi}{\prod_{J \in \mathcal J} \phi_J(\mat x_J)^{\theta_j}}  \Bigg\}
% & \Big[ \parbox{1.5in}{\centering%
% cpds $\bp$ correspond\\ to factors $\phi_j$}	\Big]\\
&\Big[ \text{definition of $\Pr_\Psi$} \Big]\\
\Ex_{\mat x \sim \mu}\! \Bigg\{ \sum_J \bigg[\,
\!\theta_J \log \frac{1}{\phi_J(\mat x_J)}
\bigg] - \sum_{J \in \mathcal J} \left[\theta_J \log \frac{1}{\phi_J(\mat x_J)} \right]
- \log Z_\Psi \Bigg\} \\
&= \Ex_{\mat x \sim \mu} [- \log Z_\Psi] \\
&= - \log Z_\Psi & \Big[~\text{$Z_\Psi$ is constant in $\mat x$}~\Big]
\end{align*}

Since $p$ has high confidence, and ${\tt T}$ is always equal to ${\tt t}$, the only joint distribution on $(X,{\tt T})$ with finite score is $\mu(X, {\tt T}) = p(X) \delta_{{\tt t}}({\tt T})$. We compute its score directly:

\begin{align*}
\aar*{\!\begin{tikzpicture}[center base]
\node[dpad0] (X) at (0,0) {$X$};
\node[dpad0] (2) at (1.1,0) {$\Truth$};
\draw[arr2] (X) to
node[above, pos=0.4,inner sep=2pt]{$\hat c$}
% node[below, pos=0.4, inner sep=2pt]{${\color{gray}\scriptstyle(\beta)}$}
\draw[arr2, <-] (X) to
node[above, pos=0.6, inner sep=2pt]{$p$}
node[below, pos=0.6, inner sep=2pt]
+(-1, 0);
\draw[arr2, <<-] (2) to
node[above, inner sep=2pt, pos=0.6]
\end{tikzpicture}\!}
%oli3: added missing \log below
&= \Ex_{\mu} \log \frac{\mu(X,{\tt T})}{\hat c({\tt t}\,|X)}
= \Ex_{p} \log \frac{1}{\hat c({\tt t}\,|X)}
= \Ex_{p} \log \frac{1}{\exp(-c(X))} \\
% \beta \Ex_{x\sim p}c(x).
&= \Ex_p \log \exp(c(X))
= \Ex_p c(X) =
\Ex_{x\sim p}\! c(x).
\end{align*}

§.§ Additional Proofs for Unnumbered Claims

§.§.§ Details on the Data Processing Inequality Proof

ci2/.style=inner sep=2pt, align=center

pstyle/.style=line width=0.9pt, pcolor!!black
qstyle/.style=line width=1.3pt, qcolor!!black
pqstyle/.style=line width=1.5pt,pcolor!50!qcolor!!black

We now provide more details on the proof of the Data Processing Equality that appeared in <Ref> of the main text. We repeat it now for convenience, with labeled PDGs ($\dg M_1, \ldots, \dg M_5$) and numbered (in)equalities.

% - \log \Pr\nolimits_{p,d}(X\!=\!x) ~=\qquad&\\
\aar*{\!\begin{tikzpicture}[center base]
\node[dpad0] (X) {$X$};
\draw[arr2, <-,qstyle] (X) --
% node[above,pos=0.6]{$q^{{\color{gray}(s)}}$}
node[below, pos=0.65,ci2] {\qlabel}
++(1.1, 0);
\draw[arr2, <-,pstyle] (X) --
% node[above,pos=0.6]{$p^{{\color{gray}(r)}}$}
node[below, pos=0.65, ci2] {\plabel}
++(-1.1, 0);%
\end{tikzpicture}\!}
\aar**{\!\!\begin{tikzpicture}[center base]
\node[dpad0] (X) {$X$};
\node[dpad0,above=.8 of X,align=center] (Y) {$Y$};
\draw[arr2, <-,qstyle] (X) --
node[below, pos=0.65,ci2] {\qlabel}
++(1.1, 0);
\draw[arr2, <-,pstyle] (X) --
node[below, pos=0.65,ci2] {\plabel}
++(-1.1, 0);%
\draw[arr2, pqstyle] (X) --
node[left,pos=0.45,inner sep=1pt]{$f$}
node[right, pos=0.45, inner sep=1.5pt, align=center] % below,rotate=90
\end{tikzpicture}\!\!}
\aar**{\!\!\begin{tikzpicture}[center base]
\node[dpad0] (X1) {$X_1$};
\node[dpad0, right=0.6 of X1] (X2) {$X_2$};
\node[dpad0,above=.8 of {$(X1)!.5!(X2)$},align=center] (Y) {$Y$};
\draw[arr2, -, double equal sign distance] (X1) to (X2);
\draw[arr2, <-,qstyle] (X2) --
node[below, pos=0.65,ci2] {\qlabel}
++(1.1, 0);
\draw[arr2, <-,pstyle] (X1) --
node[below, pos=0.65,ci2] {\plabel}
++(-1.1, 0);%
\draw[arr2,pstyle] (X1) to[bend left=40]
node[above left, pos=0.35, inner sep=1pt]{$f$}
node[below right=0 and 0, pos=0.45, inner sep=0pt, align=center] {\plabel}
\draw[arr2,qstyle] (X2) to[bend right=40]
node[above right, pos=0.35, inner sep=1pt]{$f$}
node[below left=0 and 0, pos=0.45, inner sep=0pt, align=center] {\qlabel}
\end{tikzpicture}\!\!}
\aar**{\!\!\begin{tikzpicture}[center base]
\node[dpad0] (X1) {$X_1$};
\node[dpad0, right=0.65 of X1] (X2) {$X_2$};
\node[dpad0,above=.75 of {$(X1)!.5!(X2)$},align=center] (Y) {$Y$};
\draw[arr2, <-,qstyle] (X2) --
node[below, pos=0.65,ci2] {\qlabel}
++(1.1, 0);
\draw[arr2, <-,pstyle] (X1) --
node[below, pos=0.65,ci2] {\plabel}
++(-1.1, 0);%
\draw[arr2,pstyle] (X1) to[bend left=30]
node[above left, pos=0.35, inner sep=1pt]{$f$}
node[below right=0 and 0, pos=0.45, inner sep=0pt, align=center] {\plabel}
\draw[arr2,qstyle] (X2) to[bend right=30]
node[above right, pos=0.35, inner sep=1pt]{$f$}
node[below left=0 and 0, pos=0.45, inner sep=0pt, align=center] {\qlabel}
\end{tikzpicture}\!\!}
\aar*{\!\begin{tikzpicture}[center base]
\node[dpad0] (X) {$X$};
\draw[arr2, <-,qstyle] (X) --
node[above,pos=0.7,ci2]{$ f\!\circ\! q$}
node[below, pos=0.65,ci2] {\qlabel}
++(1.1, 0);
\draw[arr2, <-,pstyle] (X) --
node[above,pos=0.6,ci2]{$ f\!\circ\! p$}
node[below, pos=0.65,ci2] {\plabel}
++(-1.1, 0);%
\end{tikzpicture}\!}
% \end{equation*}
% \]

$\dg M_1$        $\dg M_2$        $\dg M_3$        $\dg M_4$        $\dg M_5$

We now enumerate the (in)equalities to prove them.

* Let $\mu(X)$ denote the (unique) optimal distribution for $\dg M_1$.

Now, the joint distribution $\mu(X,Y) := \mu(X) f(Y|X)$ has incompatibility with $\dg M_2$ equal to

\begin{align*}
\Inc_{\dg M_2}(\mu(X,Y)) &=  \beta \kldiv{\mu(X)}{p(X)} +
\zeta \kldiv{\mu(X)}{q(X)} + (\beta\!+\!\zeta)\Ex_{x\sim\mu}\big[ \kldiv{\mu(Y|x)}{f(Y|x)} \big] \\
&= \Inc_{\dg M_1}(\mu(X)) + (\beta\!+\!\zeta) \Ex_{x \sim \mu}\kldiv{\mu(Y|x)}{f(Y|x)}
% & \hspace{-2in}{\color{gray}\Big[\text{as $\mu(X)$ is optimal for $\dg M_1$}\Big]}
\\
&= \aar{\dg M_1}
& \hspace{-2in}{\color{gray}\Big[\begin{array}{c}
\text{as $\mu(Y|x) = f(Y|x)$ wherever $\mu(x)>0$,}\\
\text{and $\mu(X)$ minimizes $\Inc_{\dg M_1}$}
\end{array}\Big]}
\end{align*}

So $\mu(X,Y)$ witnesses the fact that $\aar{\dg M_2} \le \Inc_{\dg M_2}(\mu(X,Y)) = \aar{\dg M_1}$.

Furthermore, every joint distribution $\nu(X,Y)$ must have at least this incompatibility,

as it must have some marginal $\nu(X)$, which, even by itself, already gives rise to incompatibility of magnitude $\Inc_{\dg
M_1}(\nu(X)) \ge \Inc_{\dg M_1}(\mu(X)) = \aar{\dg M_1} $.

And since this is true for all $\nu(X,Y)$, we have that $\aar{\dg M_2} \ge \aar{\dg M_1}$. So $\aar{\dg M_2} = \aar{\dg M_1}$.

* The equals sign in $\dg M_3$ may be equivalently interpreted as a cpd $\mathit{eq}_{}(X_1|X_2) := x_2 \mapsto \delta_{x_2}(X_1)$, a cpd
$\mathit{eq'}_{}(X_2|X_1) := x_1 \mapsto \delta_{x_1}(X_2)$, or both at once; in each case, the effect is that a joint distribution $\mu$ with support on an outcome for which $X_1 \ne X_2$ gets an infinite penalty, so a minimizer $\mu(X_1,X_2,Y)$ of $\Inc{\dg M_3}$ must be isomorphic to a distribution $\mu'(X,Y)$.

Furthermore, it is easy to verify that $\Inc_{\dg M_2}(\mu'(X,Y)) = \Inc_{\dg M_3}(\mu(X,X,Y))$. More formally, we have:

\begin{align*}
\aar{\dg M_3} &= \inf_{\mu(X_1,X_2, Y)} \Ex_\mu\left[
\beta \log \frac{\mu(X_1)}{p(X_1)}
+ \zeta \log \frac{\mu(X_2)}{q(X_2)}
+ \beta \log\frac{\mu(Y|X_1)}{f(Y|X_1)}
+ \zeta \log\frac{\mu(Y|X_2)}{f(Y|X_2)}
+ \log \frac{\mu(X_1|X_2)}{\mathit{eq}(X_1,X_2)}
\right]
\intertext{but if $X_1$ always equals $X_2$ (which we call simply $X$), as it must for the optimal $\mu$, this becomes}
&= \inf_{\mu(X_1=X_2=X, Y)} \Ex_\mu\left[
\beta \log \frac{\mu(X)}{p(X)}
+ \zeta \log \frac{\mu(X)}{q(X)}
+ \beta \log\frac{\mu(Y|X)}{f(Y|X)}
+ \zeta \log\frac{\mu(Y|X)}{f(Y|X)}
\right] \\
&= \inf_{\mu(X, Y)} \Ex_\mu\left[
\beta \log \frac{\mu(X)}{p(X)}
+ \zeta \log \frac{\mu(X)}{q(X)}
+ (\beta\!+\!\zeta) \log\frac{\mu(Y|X)}{f(Y|X)}
\right] \\
&= \inf_{\mu(X,Y)} \Inc_{\dg M_2}(\mu)\\
&= \aar{\dg M_2}.
% \Inc_{\dg M_2=3}(\mu'(X,Y))
\end{align*}

* Eliminating the edge or edges enforcing the equality $(X_1 = X_2)$ cannot increase inconsistency, by <Ref>.

* Although this final step of composing the edges with shared confidences looks intuitively like it should be true (and it is!), its proof may not be obvious.

We now provide a rigorous proof of this equality.

To ameliorate subscript pains, we henceforth write $X$ for $X_1$, and $Z$ for $X_2$.
We now compute:

\begin{align*}
\aar{\dg M_4}
&= \inf_{\mu(X,Z,Y)}
\Ex_\mu \left[
\beta \log \frac{\mu(X)\, \mu(Y|X)}{p(X)\,f(Y|X)}
+ \zeta \log \frac{\mu(Z)\, \mu(Y|Z)}{q(Z)\,f(Y|Z)}
\right]\\
&= \inf_{\mu(X,Z,Y)}
\Ex_\mu \left[
\beta \log \frac{\mu(Y)\, \mu(X|Y)}{p(X)\,f(Y|X)}
+ \zeta \log \frac{\mu(Y)\, \mu(Z|Y)}{q(Z)\,f(Y|Z)}
\right] & \text{[apply Bayes Rule in numerators]}
\end{align*}

By the chain rule, every distribution $\mu(X,Z,Y)$ may be specified as $\mu(Y)\mu(X|Y)\mu(Z|X,Y)$, so we can rewrite the formula above as

\begin{equation*}
\aar{\dg M_4}
\inf_{\mu(Y)} \inf_{\mu(X|Y)} \inf_{\mu(Z|Y,X)}
\Ex_{y \sim \mu(Y)} \Ex_{x \sim \mu(X|y)} \Ex_{z \sim \mu(Z|y,x)} \left[
\beta \log \frac{\mu(y)\, \mu(x\,|\,y)}{p(x)\,f(y\,|\,x)}
+ \zeta \log \frac{\mu(y)\, \mu(z\,|\,y)}{q(z)\,f(y\,|\,z)}
\right], % \\
\end{equation*}

where $\mu(Z|Y)$ is the defined in terms of the primitives $\mu(X|Y)$ and $\mu(Z|X,Y)$ as $\mu(Z|Y) := y\mapsto \Ex_{x\sim \mu(X|y)} \mu(Z|y,x)$, and is a valid cpd, since it is a mixture distribution.

Since the first term (with $\beta$) does not depend on $z$, we can take it out of the expectation, so

\begin{align*}
\aar{\dg M_4}
&= \inf_{\mu(Y)} \inf_{\mu(X|Y)} \inf_{\mu(Z|Y,X)}
\Ex_{y \sim \mu(Y)} \Ex_{x \sim \mu(X|y)} \left[
% \beta \log \frac{\mu(Y)\, \mu(X|Y)}{p(X)\,f(Y|X)}
\beta \log \frac{\mu(y)\, \mu(x\,|\,y)}{p(x)\,f(y\,|\,x)}
+~~ \zeta~ \Ex_{\substack{\vphantom{|}\\\mathclap{z\sim\mu(Z|y,x)}}}
% + \zeta \!\!\!\Ex_{\substack{\vphantom{|}\\{z\sim\mu(Z|y,x)}}} \!\!\!
\Big[ \log \frac{\mu(y)\, \mu(z\,|\,y)}{q(z)\,f(y\,|\,z)} \Big]
\right]; \\
\intertext{we can split up $\Ex_{\mu(X|y)}$ by linearity of expectation, to get}
\aar{\dg M_4}
&= \inf_{\mu(Y)} \inf_{\mu(X|Y)} \inf_{\mu(Z|Y,X)}
\Ex_{y \sim \mu(Y)} \left[
% \beta \log \frac{\mu(Y)\, \mu(X|Y)}{p(X)\,f(Y|X)}
\beta \!\! \Ex_{\substack{\vphantom{x}\\x\sim\mu(X|y)}}\!\!\Big[
\log \frac{\mu(y)\, \mu(x\,|\,y)}{p(x)\,f(y\,|\,x)} \Big]
+ \zeta\!\! \Ex_{\substack{x \sim \mu(X|y)\\ z\sim\mu(Z|y,x)}}\!\!\Big[
\log \frac{\mu(y)\, \mu(z\,|\,y)}{q(z)\,f(y\,|\,z)} \Big]
\right]
% &\text{[linearity of expectation]}
% \\
\end{align*}

Note that the quantity inside the second expectation does not depend on $x$. Therefore,
the second expectation is just an explicit way of sampling $z$ from the mixture
distribution $\Ex_{x \sim \mu(X|y)} \mu(Z|x,y)$, which is the definition of $\mu(Z|y)$.
Once we make this replacement, it becomes clear that the only feature of $\mu(Z|Y,X)$ that
matters is the mixture $\mu(Z|Y)$. Simplifying the second expectation in this way, and replacing the infemum over $\mu(Z|X,Y)$ with one over $\mu(Z|Y)$ yields:

\begin{equation*}
\aar{\dg M_4}
= \inf_{\mu(Y)} \inf_{\mu(X|Y)} \inf_{\mu(Z|Y)}
\Ex_{y \sim \mu(Y)} \left[
% \beta \log \frac{\mu(Y)\, \mu(X|Y)}{p(X)\,f(Y|X)}
\beta \!\! \Ex_{\substack{\vphantom{x}\\x\sim\mu(X|y)}}\!\!\Big[
\log \frac{\mu(y)\, \mu(x\,|\,y)}{p(x)\,f(y\,|\,x)} \Big]
+ \zeta\!\! \Ex_{\substack{\vphantom{|}\\ z\sim\mu(Z|y)}}\!\!\Big[
\log \frac{\mu(y)\, \mu(z\,|\,y)}{q(z)\,f(y\,|\,z)} \Big]
\right]
\end{equation*}

Now, a cpd $\mu(X|Y)$ is

[modulo measurability concerns that do not affect the infemum; see <Ref>]
a (possibly different) distribution $\nu_y(X)$ for every value of $Y$.
Observe that, inside the expectation over $\mu(Y)$, the cpds $\mu(X|Y)$ and $\mu(Z|Y)$ are used only for the present value of $y$, and do not reference, say,  $\mu(X|y')$ for $y'\ne y$.
Because there is no interaction between the choice of cpd $\mu(X|y)$ and $\mu(X|y')$, it is not necessary to jointly optimize over entire cpds $\mu(X|Y)$ all at once.
Rather, it is equivalent to to take the infemum over $\nu(X)$, separately for each $y$.
Symmetrically, we may as well take the infemum over $\lambda(Z)$ separately for each $y$, rather than jointly finding the optimal $\mu(Z|Y)$ all at once.
Operationallly, this means we can pull the infema inside the expectation over $Y$. And since the first term doesn't depend on $Z$ and the second doesn't depend on $X$, we get:

\begin{equation*}
\aar{\dg M_4}
= \inf_{\mu(Y)}
\Ex_{y \sim \mu(Y)} \left[
% \beta \log \frac{\mu(Y)\, \mu(X|Y)}{p(X)\,f(Y|X)}
\inf_{\nu(X)} \beta  \Ex_{\nu(X)} \Big[
\log \frac{\mu(y)\, \nu(X)}{p(X)\,f(y\,|X)} \Big]
+ \inf_{\lambda(Z)} \zeta  \Ex_{\lambda(Z)} \Big[
\log \frac{\mu(y)\, \lambda(Z)}{q(Z)\,f(y\,|Z)} \Big]
\right]
\end{equation*}

Next, we pull the same trick we've used over and over: find constants so that we can regard the dependence as a relative entropy with respect to the quantity being optimized.

Grouping the quantities apart from $\nu(X)$ on the left term and normalizing them (and analogously for $\lambda(Z)$ on the right), we find that

\begin{equation*}
\aar{\dg M_4}
= \inf_{\mu(Y)}
\Ex_{y \sim \mu(Y)} \left[
\begin{array}{l}
\beta \inf_{\nu(X)}
\kldiv* {\nu(X)}{\frac{1}{C_1(y)} p(X)\frac{f(y|X)}{\mu(y)}}
- \beta \log C_1(y) \\
+ \zeta \inf_{\lambda(Z)}
\kldiv* {\lambda(Z)}{\frac{1}{C_2(y)} q(Z)\frac{f(y|Z)}{\mu(y)}}
- \zeta \log C_2(y)
\end{array}
\right],
\end{equation*}

\[
C_1(y) = \sum_x p(x)\frac{f(y|x)}{\mu(y)} = \frac{1}{\mu(y)} \Ex_{p(X)} f(y|X)
\qquad\text{and}\qquad
C_2(y) = \sum_z q(z)\frac{f(y|z)}{\mu(y)} = \frac{1}{\mu(y)} \Ex_{q(Z)} f(y|Z)
\]

are the constants required to normalize the distributions. Both relative entropies are minimized when their arguments match, at which point they contribute zero, so we have

\begin{align*}
\aar{\dg M_4}
&= \inf_{\mu(Y)}
\Ex_{y \sim \mu(Y)} \left[
% \beta \inf_{\nu(X)}
% \kldiv* {\nu(X)}{\frac{1}{C_1(y)} p(X)\frac{f(y|X)}{\mu(y)}}
\beta \log \frac1{C_1(y)}
% + \zeta \inf_{\lambda(Z)}
% \kldiv* {\lambda(Z)}{\frac{1}{C_2(y)} q(Z)\frac{f(y|Z)}{\mu(y)}}
+ \zeta \log \frac1{C_2(y)}
\right]\\
&= \inf_{\mu(Y)}
\Ex_{y \sim \mu(Y)} \left[
\beta \log \frac{\mu(y)}{\Ex_{p(X)} f(y|X)}
+ \zeta \log \frac{\mu(y)}{\Ex_{q(Z)} f(y|Z)} \right] \\
&= \inf_{\mu(Y)} \Ex_{\mu} \Big[ \beta \kldiv{\mu}{f \circ p} + \zeta \kldiv{\mu}{f\circ q}  \Big] \\
&= \aar{\dg M_5}.
\end{align*}

§.§.§ Details for Claims made in Section 8

First, the fact that

\[
% \dg M_1 := ~
\mathcal L_1 = \lambda_\dsymb\mathcal L_\datsymb + \lambda_\ssymb \mathcal L_\simsymb = \aar**{
\begin{tikzpicture}[center base]
% \node[dpad0] (Z) at (-0.2,0) {$Z$};
% \node[tpt={z0|$0$}] at (-0.5,0.1) {};
% \node[tpt={z1|$1$},right=0.15 of z0]{};
\node[tpt={z0|\simsymb}] at (-0.5,0.1) {};
\node[tpt={z1|\datsymb},right=0.35 of z0]{};
\node[Dom={$Z$[label distance=-2.5ex, xshift=1.0em] (Z)
around {\lab{z0}\lab{z1}}},yshift=0.2em ] {};

% \node[dpad0,align=center] (XY) at (1.8,0) {$XY$}; %{$X$\\[-0.3ex]$Y$};
\node[dpad0] (X) at (2.4, 0.6) {$X$};
\node[dpad0] (Y) at (2.4, -0.6) {$Y$};
\coordinate (xyz) at (1.9, 0);
\draw[arr1, <-] (Z) to
% node[above, pos=0.6]{$\hat\lambda$}
node[above, pos=0.6]{$\lambda$}
% node[above, pos=0.6]{$\frac{1}{\lambda_0+\lambda_1}[\lambda_0, \lambda_1]$}
node[below,inner sep=1pt, pos=0.6]{${\color{gray}\scriptstyle( \infty )}$}
+(-1.5, 0);
% \node at (-1,-0.6) {\small where $\lambda(Z) = \frac{\lambda_Z}{\lambda_0+\lambda_1}$};
% \node at (0,-0.6) {\small where $\lambda(Z) \propto \lambda_Z$};
\draw[arr1] (X) to node[right,pos=0.4]{$h$} (Y);
\draw[arr,-,shorten >=0pt] (Z) to[bend left=0, shorten >=0pt]
% node[fill=white, inner sep=0pt, pos=0.55]
% node[inner sep=1pt, pos=0.55]
node[above, inner sep=1pt, pos=0.55]
% {$Z?d:s$}
% {$\begin{bmatrix}d \text{ if } Z\\[-0.3ex] s \text{ else}\end{bmatrix}$}
% {$\begin{matrix}d \text{ if } Z\!\!=\!\!1\\[-0.3ex]
% 	% s \text{ else}\end{matrix}$}
% 	\text{else }s\end{matrix}$}
{$\begin{matrix}\datsymb \mapsto d \\[-0.6ex]
\simsymb \mapsto s \end{matrix}$}
% node[above, inner sep=2pt, pos=0.68]
% 	{${\color{gray}\scriptscriptstyle(r)}$}
node[below,inner sep=1pt]{${\color{gray}\scriptstyle( \infty )}$}
\draw[arr2, shorten <=0pt] (xyz) to (X);
\draw[arr2, shorten <=0pt] (xyz) to (Y);
\end{tikzpicture}}
% (\lambda_\ssymb+\lambda_\dsymb),
\]

where $\lambda(Z=\simsymb) = \lambda_\ssymb$ and
$\lambda(Z=\datsymb) = \lambda_\dsymb$
is immediate.
The two cpds with infinite confidence ensure that the only joint distribution with a finite score is $\lambda_\ssymb s + \lambda_\dsymb d$, and the inconsistency with $h$ is its surprisal, so the inconsistency of this PDG is

\begin{align*}
\Ex_{\lambda_\ssymb s + \lambda_\dsymb d} \Big[\log \frac{1}{h(Y|X)}\Big]
= - \lambda_\ssymb \Ex_{s} [\log {h(Y|X)}] - \lambda_\dsymb  \Ex{d}[\log h(Y|X)]
= \lambda_\dsymb\mathcal L_\datsymb + \lambda_\ssymb \mathcal L_\simsymb
= \mathcal L_1,
\quad\text{as promised.}
\end{align*}

The second correspondence is the least straightforward. Let  $C = \int sd$ be the normalization constant required to normalize the joint density $sd$. We claim that, for large fixed $\gamma$, we have

\[
% \dg M_2 :=
\mathcal L_2 \approx
\aar**{
\begin{tikzpicture}[center base]
\node[dpad0] (X) at (0, 0.6) {$X$};
\node[dpad0] (Y) at (0, -0.6) {$Y$};
\draw[arr1] (X) to node[left, pos=0.4, inner sep=1pt]{$h$}
% node[below=0pt,inner sep=1pt,rotate=90]{${\color{gray}\scriptstyle(\!\alpha{:}0\!)}$}
%oli: NO NEED!
% node[right=0pt,inner sep=1pt]{${\color{gray}\scriptstyle
% 	\renewcommand{\arraystretch}{.7}
% 	\big(\begin{matrix}
% 		\scriptstyle \alpha : 0 \\ \scriptstyle \beta : 1
% 	\end{matrix}
% 	\big)}$}

\coordinate (d0) at (1.8, 0);
\coordinate (dmid) at (0.9, 0);
\coordinate (s0) at (-1.8, 0);
\coordinate (smid) at (-0.9, 0);

\draw[arr,->,shorten <=0pt] (dmid) to[bend right=25] (X);
\draw[arr,->,shorten <=0pt] (dmid) to[bend left=25] (Y);
\draw[arr1,-,shorten <=0pt] (dmid) to
node[below, inner sep=2pt]{${\color{gray}\scriptstyle
\renewcommand{\arraystretch}{.7}
\big(\begin{matrix}
\scriptstyle\alpha: 1 \\[-0.2ex] \scriptstyle\beta: \gamma
\end{matrix} \big)}$}
node[above] {$d$}
\draw[arr,->,shorten <=0pt] (smid) to[bend left=25] (X);
\draw[arr,->,shorten <=0pt] (smid) to[bend right=25] (Y);
\draw[arr1,-,shorten <=0pt] (smid) to
node[below, inner sep=2pt]{${\color{gray}\scriptstyle
\renewcommand{\arraystretch}{.7}
\big( \begin{matrix}
\scriptstyle \alpha: 1 \\[-0.2ex] \scriptstyle \beta: \gamma
\end{matrix} \big)}$}
\end{tikzpicture}}\Bigg._{\!\!\!\gamma}
% - k \log Z_{sd} + H
% - k \log C,
+ \mathit{const},
% \overbrace{ - k \log C, }^{\text{normalization constant for $sd$}}
\]

where $\mathit{const}$ does not depend on $h$. To see this, let $\dg M_2$ be the PDG above, and compute

\begin{align*}
% \aar**{\!\!\!
% \begin{tikzpicture}[center base]
% 	\node[dpad0] (X) at (0, 0.6) {$X$};
% 	\node[dpad0] (Y) at (0, -0.6) {$Y$};
% 	\draw[arr1] (X) to node[left, pos=0.4, inner sep=1pt]{$h$} (Y);
% 	\coordinate (d0) at (1.8, 0);
% 	\coordinate (dmid) at (0.9, 0);
% 	\coordinate (s0) at (-1.8, 0);
% 	\coordinate (smid) at (-0.9, 0);
% 	\draw[arr,->,shorten <=0pt] (dmid) to[bend right=25] (X);
% 	\draw[arr,->,shorten <=0pt] (dmid) to[bend left=25] (Y);
% 	\draw[arr1,-,shorten <=0pt] (dmid) to
% 		node[below, inner sep=2pt]{${\color{gray}\scriptstyle
% 			\renewcommand{\arraystretch}{.7}
% 			\big(\begin{matrix}
% 				\scriptstyle\alpha: 1 \\[-0.2ex] \scriptstyle\beta: \gamma
% 			\end{matrix} \big)}$}
% 		node[above] {$d$}
% 		(d0);
% 	%
% 	\draw[arr,->,shorten <=0pt] (smid) to[bend left=25] (X);
% 	\draw[arr,->,shorten <=0pt] (smid) to[bend right=25] (Y);
% 	\draw[arr1,-,shorten <=0pt] (smid) to
% 		node[below, inner sep=2pt]{${\color{gray}\scriptstyle
% 			\renewcommand{\arraystretch}{.7}
% 			\big( \begin{matrix}
% 				\scriptstyle \alpha: 1 \\[-0.2ex] \scriptstyle \beta: \gamma
% 			\end{matrix} \big)}$}
% 		node[above]{$s$}
% 		(s0);
% \end{tikzpicture}\!\!}\Bigg._{\!\!\!\gamma}
\aar{\dg M_2}_\gamma
&= \inf_{\mu(X,Y)} \Ex_{\mu} \bigg[
\overbracket{
% \gamma \log \frac{\mu(XY)}{s(XY)}
% + \gamma \log \frac{\mu(XY)}{d(XY)}
\gamma \log \frac{\mu(XY)}{s(XY)}
\frac{\mu(XY)}{d(XY)}
+ \log \frac{\mu(Y|X)}{h(Y|X)} }^{\Inc(\mu)}
\overbracket{
%   \gamma \log \frac{1}{s(XY)}
% + \gamma \log \frac{1}{d(XY)}
\gamma \log \frac{1}{s(XY)}
\frac{1}{d(XY)}
- \gamma \log \frac{1}{\mu(XY)}
\bigg] \\
&= \inf_{\mu(X,Y)} \Ex_{\mu} \bigg[
\gamma \log \frac{\mu(XY)}{s(XY)}
\frac{\mu(XY)}{d(XY)}
\frac{1}{\mu(XY)}
\frac{1}{\mu(XY)}
\frac{\mu(XY)}{1}
+ \log \frac{\mu(Y|X)}{h(Y|X)}
\bigg] \\
&= \inf_{\mu(X,Y)} \Ex_{\mu} \bigg[
\gamma \log
\frac{\mu(XY)}{s(XY)d(XY)}
+ \log \frac{\mu(Y|X)}{h(Y|X)}
\bigg] \\
&= \inf_{\mu(X,Y)} \Ex_{\mu} \bigg[
\gamma \log
\frac{\mu(XY) C}{s(XY)d(XY)} - \gamma \log C
+ \log \frac{\mu(Y|X)}{h(Y|X)}
\bigg] \\
&= \inf_{\mu(X,Y)}
\gamma \kldiv*{\mu}{\frac1C{sd}} +
\Ex_{\mu} \bigg[ \log \frac{\mu(Y|X)}{h(Y|X)}\bigg]
- \gamma \log C \\
\end{align*}

$\thickD$ is $(\gamma m)$-strongly convex in a region around its minimizer for some $m>0$ that depends only on $s$ and $d$.
Together with our assumption that $h$ is positive, we find that when $\gamma$ becomes large,

the first term dominates, and the optimizing $\mu$ quickly approaches the normalized density $\nu := \frac1Csd$.
Plugging in $\nu$, we find that the value of the infemum approaches

\begin{align*}
\aar{\dg M_2} &\approx \Ex_{\nu} \bigg[ \log \frac1{h(Y|X)} \bigg] - H_{\nu}(Y|X) - \gamma \log C \\
&= \int_{XY} \frac1C \log \frac1{h(Y|X)} s(X,Y) d(X,Y) \quad- H_{\nu}(Y|X) - \gamma \log C \\
&= \frac{1}{C} \Ex_{s} \bigg[ d(X,Y) \log \frac1{h(Y|X)} \bigg] - H_{\nu}(Y|X) - \gamma \log C  \\
&= \frac{1}{C} \mathcal L_2 - H_{\nu}(Y|X) - \gamma \log C,   \\[1ex]
% \implies\qquad
\text{and therefore}\qquad
\mathcal L_2 &= C \aar{\dg M_2} + C \H_\nu(Y|X) - \gamma \, C \log C
\\ &= C \aar{\dg M_2} + \mathit{const}.
\end{align*}

Finally, we turn to

\[
% \dg M_3 :=
\mathcal L_3 := \aar**{
\begin{tikzpicture}[center base]
\node[dpad0] (X) at (0, 0.6) {$X$};
\node[dpad0] (Y) at (0, -0.6) {$Y$};
\draw[arr1] (X) to node[left=0pt,pos=0.4, inner sep=1pt]{$h$} (Y);

% \coordinate (d0) at (1.3, 0);
% \coordinate (s0) at (-1.3, 0);
% \node[above left=1pt and 0.5em of d0] {$d$};
% \node[below left=0pt and 0.2em of d0]{${\color{gray}\scriptstyle( \lambda_1 )}$};
% \node[above right=1pt and 0.5em of s0] {$s$};
% \node[below right=0pt and 0.2em of s0]{${\color{gray}\scriptstyle( \lambda_0 )}$};

\coordinate (d0) at (1.8, 0);
\coordinate (dmid) at (0.9, 0);
\coordinate (s0) at (-1.8, 0);
\coordinate (smid) at (-0.9, 0);

\draw[arr,->,shorten <=0pt] (dmid) to[bend right=25] (X);
\draw[arr,->,shorten <=0pt] (dmid) to[bend left=25] (Y);
\draw[arr1,-,shorten <=0pt] (dmid) to
node[below, inner sep=2pt]{${\color{gray}\scriptstyle(\lambda_{\dsymb})}$}
node[above] {$d$}
\draw[arr,->,shorten <=0pt] (smid) to[bend left=25] (X);
\draw[arr,->,shorten <=0pt] (smid) to[bend right=25] (Y);
\draw[arr1,-,shorten <=0pt] (smid) to
node[below, inner sep=2pt]{${\color{gray}\scriptstyle(\lambda_{\ssymb})}$}
% \unmergearr{s0}XY
% \unmergearr{d0}XY
% 	\draw[arr,-,shorten >=0pt] (Z) to[bend left=0, shorten >=0pt]
% 		node[fill=white, inner sep=0pt, pos=0.55]
% 			% {$Z?d:s$}
% 			{$\begin{bmatrix}d\\[-0.3ex] s\end{bmatrix}$}
% 		% node[above, inner sep=2pt, pos=0.68]
% 		% 	{${\color{gray}\scriptscriptstyle(r)}$}
% 		(xyz);
% 	\draw[arr2, shorten <=0pt] (xyz) to (X);
% 	\draw[arr2, shorten <=0pt] (xyz) to (Y);
% ode[dpad0] (X) {};
\end{tikzpicture}}.
\]

To see the why the optimal distribution $\mu^*(XY)$ is the $\lambda$-weighted geometric mean of $s$ and $d$

, let us first consider the same PDG, except without $h$.
From <Ref>, we have this loss without $h$ in closed form, and from the proof of <Ref>, we see that the optimizing distribution in this case is
the $\lambda$-weighted geometric distribution $\mu^* \propto s(XY)^{\lambda_\ssymb} d(XY)^{\lambda_\dsymb}$.
Now (<Ref>), including $h$ cannot make the PDG any less inconsistent. In particular, by choosing

\[
h^*(Y|X) := \mu^*(Y|X) \propto (Y|X)^{\lambda_\ssymb} d(Y|X)^{\lambda_\dsymb},
\]

to be already compatible with this joint distribution, the inconsistency does not change, while choosing a different $h$ would cause the inconsistency to increase. Thus, the optimal classifier $h^*$ by this metric is indeed as we claim.  Finally, it is easy to see that this loss is calibrated: if $s = d$, then the optimal joint distribution is equal to $s$ and to $d$, and the optimal classifier is $h(Y|X) = s(Y|X) = d(Y|X)$. So $\mathcal L_3$ is calibrated.

§.§.§ Details for Claims made in Section 9

Distortion Due to Inconsistency.
In the footnote on fn:logEexp, we claimed that if the model confidence $\beta_p$ were 1 rather than $\infty$, we would have obtained an incconsistency of

$ - \log \Ex_{x\sim p} \exp(- c(x)) $,
and that the optimal distribution would not have been $p(X)$.

\begin{align*}
\aar*{\!\begin{tikzpicture}[center base]
\node[dpad0] (X) at (0,0) {$X$};
\node[dpad0] (2) at (1.1,0) {$\Truth$};
\draw[arr2] (X) to
node[above, pos=0.4,inner sep=2pt]{$\hat c$}
% node[below, pos=0.4, inner sep=2pt]{${\color{gray}\scriptstyle(\beta)}$}
\draw[arr2, <-] (X) to
node[above, pos=0.6, inner sep=2pt]{$p$}
% node[below, pos=0.6, inner sep=2pt]
% {${\color{gray}\scriptscriptstyle(\mskip-2mu\infty\mskip-2mu)}$}
+(-1, 0);
\draw[arr2, <<-] (2) to
node[above, inner sep=2pt, pos=0.6]
\end{tikzpicture}\!}
&=  \inf_{\mu(X)} \Ex_{x \sim \mu} \left[ \log \frac{\mu(x)}{p(x)}
+ \log \frac{\mu(\trut\,|\,x)}{\hat c(\trut\,|\,x)} \right]\\
&=  \inf_{\mu(X)} \Ex_{x \sim \mu} \left[ \log \frac{\mu(x)}{p(x)}
+ \log \frac{1}{\hat c(\trut\,|\,x)} \right] \\
&=  \inf_{\mu(X)} \Ex_{x \sim \mu} \left[ \log \frac{\mu(x)}
{p(x) \exp(-c(x))}\cdot\frac{Z}{Z} \right] \\
\intertext{\raggedleft where $Z = \sum_x p(x) \exp(-c(x)) = \Ex_p \exp(-c(X))$ is the constant required to normalize the distribution
% $\lambda(X) := \frac{1}{Z}p(X)\exp(-c(X)).
&=  \inf_{\mu(X)} \kldiv*{\mu}{\frac{1}{Z}\,p(X)\exp(-c(X))} - \log Z \\
&=  - \log Z \\
&=  - \log \Ex_{x\sim p} \exp(-c(x))
\end{align*}

as promised.
Note also that in the proof, we showed that the optimal distribution is proportional to $p(X) \exp(-c(X))$ which means that it equals $p(X)$ if and only if $c(X)$ is constant in $X$.

Enforcing the Qualitative Picture.
We also claimed without careful proof in <Ref> that, if $\alpha_h = \alpha_{\datadist\xysamp} = 1$, then

\begin{equation*}
\lim_{\gamma\to\infty}
\left.
\aar**{\begin{tikzpicture}[center base]
\begin{scope}[xscale=1.2]
\node[dpad0] (X) at (0.3,0) {$X$};
\node[dpad0] (Yt) at (1,1) {$Y$};
\node[dpad0,align=center] (Yp) at (1.4,0) {$\vphantom{Y}\smash{Y'}$};
\node[dpad0] (2) at (2,1) {$\Truth$};
\coordinate (dstart) at (-0.1,0.9);
\end{scope}
\unmergearr[arr1]{dstart}{X}{Yt}
\node[above=2pt of center-dstartXYt, xshift=-2pt] {$\datadist\xysamp$};
\node[below right=2.0pt and -0.4pt of center-dstartXYt, inner sep=0pt, rotate=25]
\mergearr[arr2]{Yt}{Yp}{2}
\node[above=2pt of center-YtYp2] {$\hat\ell$};
\draw[arr2] (X) to
node[above, inner sep=2pt,pos=0.4] {$h$}
node[below, inner sep=2pt,pos=0.4]
\draw[arr2, <<-] (2) to
node[right, inner sep=2pt, pos=0.6]
% +(0,1);
\end{tikzpicture}}\right._{\!\!\!\gamma}
= \quad\;\;\mathop{\scalebox{1.2}{$\Ex$}}\limits_{\substack{%
\vphantom{x}\\
\mathllap{(x,y)} \sim \mathrlap{\datadist\xysamp} \\
\mathllap{y'} \sim \mathrlap{p(Y'|\,x)}} }
\;\big[\ell(y,y')\big]
\end{equation*}

Why is this? For such a setting of $\alpha$, which intuitively articulates a causal picture where $X,Y$ is generated from $\datadist\xysamp$, and $Y'$ generated by $h(Y'|X)$, the information deficiency $\IDef{\dg S}(\mu(X,Y,Y'))$ of a distribution $\mu$ is

\begin{align*}
\IDef{\dg S}(\mu(X,Y,Y')) &= -\H_\mu(X,Y,Y') + \H(X,Y) + \H(Y'|X) \\
&= \H_\mu(Y'|X) - \H_\mu(Y' | X, Y)  \\
&= \I_\mu(Y;Y'|X).
\end{align*}

Both equalities of the derivation above standard information theoretic identities [See, for instance,][]mackay2003information, and the final quantity $\I_{\mu}(Y;Y'|X)$ is the conditional mutual information between $Y$ and $Y'$ given $X$, and is a non-negative number that equals zero if and only if $Y$ and $Y'$ are conditionally independent given $X$.

As a result, as $\gamma\to\infty$ any distribution that for which  $Y'$ and $Y$ are not independent given $X$ will incur infinite cost. Since the confidences in $h$ and $\datadist\xysamp$ are also infinite, so will a violation of either cpd.
There is only one distribution that has both cpds and also this independence; that distribution is $\mu(X,Y,Y') := \datadist\xysamp(X,Y)h(Y'|X)$.
Now the argument of <Ref> applies: all other cpds must be matched, and the inconsistency is the expected incompatibility of $\hat l$, which equals

\[ \quad\;\;\mathop{\scalebox{1.2}{$\Ex$}}\limits_{\substack{%
\vphantom{x}\\
\mathllap{(x,y)} \sim \mathrlap{\datadist\xysamp} \\
\mathllap{y'} \sim \mathrlap{p(Y'|\,x)}} }
\; \log\frac{1}{\hat\ell(\trut\,|y,y')}
\quad\;\;\mathop{\scalebox{1.2}{$\Ex$}}\limits_{\substack{%
\vphantom{x}\\
\mathllap{(x,y)} \sim \mathrlap{\datadist\xysamp} \\
\mathllap{y'} \sim \mathrlap{p(Y'|\,x)}} }
\;  \log\frac{1}{\exp(-\ell(y,y'))}
\quad\;\;\mathop{\scalebox{1.2}{$\Ex$}}\limits_{\substack{%
\vphantom{x}\\
\mathllap{(x,y)} \sim \mathrlap{\datadist\xysamp} \\
\mathllap{y'} \sim \mathrlap{p(Y'|\,x)}} }
\;\big[ \log \exp(\ell(y,y')) \big]
\quad\;\;\mathop{\scalebox{1.2}{$\Ex$}}\limits_{\substack{%
\vphantom{x}\\
\mathllap{(x,y)} \sim \mathrlap{\datadist\xysamp} \\
\mathllap{y'} \sim \mathrlap{p(Y'|\,x)}} }
\;\big[\ell(y,y')\big]
= \mathcal L
\]

§ MORE NOTES

§.§ Maximum A Posteriori and Priors

The usual telling of the correspondence between regularizers and priors is something like the following.

Suppose you have a parameterized family of distributions

and have observed evidence $X$, but do not know the parameter $\Theta$.
The maximum-likelihood estimate of $\Theta$ is then

\[
\theta^{\mathrm{MLE}}(X) := \arg\max_{\theta\in \Theta}  \Pr(X|\theta)
= \arg\max_{\theta\in \Theta} \log \Pr(X|\theta).
\]

The logarithm is a monotonic transformation, so it does not change the argmax, but it has
nicer properties, so that function is generally used instead. (Many of the loss functions in main body of the paper are log-likelihoods also.)

In some sense, better than estimating the maximum likelihood, is to perform a Bayesian update with the new information, to get a distribution over $\Theta$.
If that's too expensive, we could simply take the estimate with the highest posterior probability, which is called the Maximum A Posteriori (MAP) estimate.
For any given $\theta$, the Bayesian reading of Bayes rule states that

\[
\text{posterior $\Pr(\Theta | X)$} = \frac
{\text{likelihood $\Pr(X|\Theta)$}\cdot\text{prior $\Pr(\Theta)$}}{\text{evidence $\Pr(X) = \sum_{\theta'} \Pr(X|\theta')\Pr(\theta')$}}.
\]

So taking a logarithm,

\[
\text{log-posterior $\log \Pr(\Theta | X)$} = \text{log-likelihood $\log \Pr(X|\Theta)$} ~+~ \text{log-prior $\log \Pr(\Theta)$}
- \text{log-evidence $\log \Pr(X)$}.
\]

The final term does not depend on $\theta$, so it is not relevant for finding the optimal $\theta$ by this metric. Swapping the signs so that we are taking a minimum rather than a maximum, the MAP estimate is then given by

\[
\theta^{\mathrm{MAP}}(X) := \arg\min_{\theta \in \Theta} \left\{ \log \frac{1}{\Pr(X|\theta)} + \log \frac1{\Pr(\theta)} \right\}.
\]

Note that if negative log likelihood (or surprisal, $-\log \Pr(X|\theta)$) was our original loss function, we have now added an arbitrary extra term, as a function of $\Theta$, to our loss function. It is in this sense that priors classically correspond to regularizers.

§.§ Surprise

A common justification for using $\I_p(x)$ as a cost for updating a probabilistic model $p(x)$ based on an observed sample $x$, is that by minimizing it, you  “maximize the probability of seeing your data”.

[this justification should not be taken too seriously  without constraints on $p$, because the optimal value of $p$  is $\delta_x$, which does not generalize.]
But this explanation applies just as well to $-p(x)$. Why include the logarithm?
There are plenty of answers to this question; among them: $\I_p$ is convex in $p$, it decomposes products into arguably simpler sums, is more numerically stable, has a well-defended physical analog in thermodynamics, and is a primative of information theory.

For those after a quick and rigorous justification (as opposed to handwaving or a thermodynamics textbook), none of these answers are entirely satisfying.
They suggest that $\I_p$ has certain nice properties, but not that it enjoys them uniquely, or that no other loss function satisfies nicer ones.
Pedagogically speaking, the situation is more straightforward for us.
Although PDG semantics themselves require non-trivial justification

, they give us in return uniform answers to many questions, starting with:
Why use the surprise $\I_p(x)$, to measure the loss of a model $p(X)$ on sample $x$? Because it is the inconsistency of simultanously believing $X = x$ and $X \sim p$.