# Magnetovolume Effect on the First-Order Metamagnetic Transition in UTe2

Atsushi Miyake1<EMAIL_ADDRESS>Masaki Gen2 Akihiko Ikeda3 Kazumasa
Miyake4 Yusei Shimizu5 Yoshiki J. Sato5 Dexin Li5 Ai Nakamura5 Yoshiya Homma5
Fuminori Honda5,6 Jacques Flouquet7
Masashi Tokunaga1 and Dai Aoki5,7 1Institute for Solid State Physics
(ISSP)1Institute for Solid State Physics (ISSP) The University of Tokyo The
University of Tokyo Kashiwa Kashiwa Chiba 277-8581 Chiba 277-8581 Japan
2Department of Advanced Materials Science Japan
2Department of Advanced Materials Science The University of Tokyo The
University of Tokyo Kashiwa Kashiwa Chiba 277-8561 Chiba 277-8561 Japan
3Department of Engineering Science Japan
3Department of Engineering Science University of Electro-Communications
University of Electro-Communications Chofu Chofu Tokyo 182-8585 Tokyo 182-8585
Japan
4Center for Advanced High Magnetic Field Science Japan
4Center for Advanced High Magnetic Field Science Osaka University Osaka
University Toyonaka Toyonaka Osaka 560-0043 Osaka 560-0043 Japan
5Institute for Materials Research Japan
5Institute for Materials Research Tohoku University Tohoku University Oarai
Oarai Ibaraki 311-1313 Ibaraki 311-1313 Japan
6Central Institute of Radioisotope Science and Safety Japan
6Central Institute of Radioisotope Science and Safety Kyushu University Kyushu
University Fukuoka Fukuoka 819-0395 819-0395 Japan
7University Grenoble Alpes Japan
7University Grenoble Alpes CEA CEA Grenoble INP Grenoble INP IRIG-PHELIQS
IRIG-PHELIQS 38000 Grenoble 38000 Grenoble France France

###### Abstract

The link between the metamagnetic transition and novel spin-triplet
superconductivity of UTe2 was discussed thermodynamically through
magnetostriction measurements in a pulsed-magnetic field. We revealed a
discontinuous magnetostriction across the metamagnetic transition at
$\mu_{0}H_{\rm m}\approx 35$ T for the applied magnetic fields along the
crystallographic $b$ axis in the orthorhombic structure. The resultanting
volume magnetostriction of $\Delta V/V\approx-5.9\times 10^{-4}$ gives the
initial pressure dependence of $H_{\rm m}$ by employing the Clausius-
Clapeyron’s equation, which agrees with previous pressure experiments.
Further, significant anisotropic magnetostriction (AMS), derived by
subtracting the averaged linear magnetostriction, was revealed. Contrary to
the weakly field-dependent AMS along the $a$ axis, those along the $b$ and $c$
axes show strong field dependences with a similar magnitude but with opposite
signs, indicating its lattice instability. The relationship between
characteristic energy scales of magnetic fields and temperatures was discussed
in terms of the Grüneisen parameters compared to the other $f$-electron
systems. The volume shrinkage in UTe2 at $H_{\rm m}$, contrary to the volume
expansion in typical heavy fermion metamagnets, pushes to invoke the link with
the valence instability related to the itinerant-localized dual nature of the
U magnetism.

Unconventional superconductivity of UTe2 (space group: $Immm$) below $T_{\rm
sc}$ = 1.6-2.0 K is a promising candidate for a spin-triplet state [1, 2, 3].
Since its discovery few years ago, a number of experimental and theoretical
works have been intensively conducted [4]. For instance, puzzling magnetic-
field ($H$) effects have been reported from the early stage. A first-order
metamagnetic transition takes place at $\mu_{0}H_{\rm m}\approx 35$ T, when
the field is applied along the $b$ axis [5, 6, 7, 8, 9] [see Fig. 1(a)].
Notable features appear in its superconducting (SC) phase diagrams, which are
significantly sensitive to the field and its applied direction. For
$H~{}||~{}b$, $T_{\rm sc}$ decreases to $\sim$15 T, followed by a gradual
increase and a sudden suppression at $H_{\rm m}$ [10, 7, 8]. More strikingly,
another SC phase is induced above $H_{\rm m}$ for a magnetic field along
specific directions near the [011] direction [7, 8]. Such
reinforcement/reentrant behaviors of SC transition are reminiscent of spin-
triplet ferromagnetic (FM) superconductors of URhGe and UCoGe [11, 12, 13, 14,
15]. However, UTe2 does not order magnetically [16, 17, 18]. The spin
fluctuations have been investigated microscopically by NMR [19] and $\mu$SR
[16]. Direct antiferromagnetic (AFM) correlations are detected using inelastic
neutron scattering experiments [20, 21], while FM correlations have not yet
been observed experimentally up to date. These AFM and FM fluctuations may
compete for the choice of SC order parameters [22, 23].

Although it is evident that the first-order metamagnetic transition influences
the SC phases, the relation between them, particularly the origin of
metamagnetism, leaves many mysteries. For $H~{}||~{}b$, a steep increase in
magnetization from $M\approx 0.5~{}\mu_{\rm B}$/U to $M\approx 1.0~{}\mu_{\rm
B}$/U at $H_{\rm m}$ [see Fig. 1(a)] accompanied by approximately five times
enhancement of the resistivity with a clear $H$-hysteresis at low temperatures
is observed [5, 6, 8, 9]. Through magnetization and specific heat
measurements, we revealed a significant effective mass enhancement on
approaching $H_{\rm m}$ [6, 24, 9]. A similar mass enhancement was also
supported by resistivity measurements [5, 8] and explained theoretically[25].
Note that above $H_{\rm m}$, the $b$ axis becomes the easy magnetization axis,
while below $H_{\rm m}$, it is the $a$ axis with a quasi-saturation to $\sim$1
$\mu_{\rm B}$/U at $\mu_{0}H\approx 20$ T [6, 9]. The singular point is that
the slope $dM/dH$ for $H~{}||~{}b$ on both side of $H_{\rm m}$ is quasi
constant, more precisely only $\sim$14% reduction across $H_{\rm m}$, although
the jump $\Delta M$ at $H_{\rm m}$ is generally regarded as a mark of drastic
change in the localized 5$f$ contribution.

Another remarkable feature in UTe2 is an intermediate valence state at ambient
pressure revealed by uranium $L_{3}$ x-ray absorption near-edge spectroscopy
(XANES) [26] and core-level photoelectron spectroscopy measurements[27]. From
the latter experiments, the dominant 5$f^{3}$ configuration is suggested to
contribute to the uranium 5$f$ state in UTe2. Moreover, the XANES results
under pressure ($p$) revealed an apparent uranium valence shift of $\sim$0.1
across $p_{c}\approx 1.4$ GPa[26], where the SC phase disappeared, instead,
magnetically ordered phases appeard. [28, 29, 30, 31, 32, 33, 1]. This valence
shift coincides with a change in the magnetic anisotropy. The easy
magnetization axis switched from the $a$ to $b$ axes above $p_{c}$ [1]. At
$p_{c}$, the $H_{\rm m}$ was also suppressed, as shown in Fig. 2 [31, 32].
Thus, the valence shift was accompanied by the collapse of the metamagnetic
transition at $H_{\rm m}$. Such a valence instability is also expected at
$H_{\rm m}$ at ambient pressure.

In this study, we measure the linear magnetostriction[35] of UTe2 along all
the principal axes for $H~{}||~{}b$ and estimate the resultant volume
magnetostriction. As expected from thermodynamic relations, a sharp step-like
volume shrinkage across $H_{\rm m}$ is clearly observed. Using the
magnetostriction data, we discuss the metamagnetism from the thermodynamic
point of view and a possible link to the valence instability at $H_{\rm m}$ at
ambient pressure.

Single crystals of UTe2 were grown using the chemical vapor transport method
[2]. The sample used here was cut from the same piece that reported the
simultaneous measurements of magnetization and magnetocaloric effect (MCE)
[9]. Magnetostriction in pulsed-magnetic fields was measured using fiber-
Bragg-grating (FBG) and the optical filter method [36]. The bare optical
fibers with FBG were glued directly on the sample parallel and perpendicular
to the field direction along the $b$ axis, as in Ref. Gen2022. This procedure
enabled us to measure the longitudinal and transverse magnetostriction
simultaneously, the first attempt in a pulsed-field. Thus, two datasets of
linear magnetostriction along the $a$ and $b$ axes, and along the $b$ and $c$
axes, were measured to obtain the volume magnetostriction. Almost identical
results along the $b$ axis for the different experimental setups confirmed the
reliability of our magnetostriction measurements. The measurements were
performed at low temperatures down to 1.4 K, where the sample was immersed in
4He liquid [38]. Pulsed-magnetic fields up to $\sim$54 T were generated using
non-destructive pulse-magnets with typical pulses durations of $\sim$36 ms,
installed at the ISSP of the University of Tokyo.

Figure 1: (Color online) Magnetic field dependence of (a)magnetization,
(b)linear and volume ($V$) magnetostriction, and (c)anisotropic linear
magnetostriction of UTe2 at $T_{\rm ini}$ = 1.4 K for $H~{}||~{}b$. The volume
magnetostriction is determined by a summation of the linear magnetostriction,
$\Delta V/V=\sum_{i=a,b,c}(\Delta L_{i}/L_{i})$. $L_{\rm avg}$ is an average
of the volume (lattice) magnetostriction determined by $\frac{\Delta V}{3V}$.
The anisotropic magnetostriction is evaluated by $\Delta L_{i}^{\rm
AMS}/L_{i}=\Delta L_{i}/L_{i}-\Delta V/(3V)$. The $M(H)$ curve is taken from
Refs. Miyake2019 and AMiyake2021.

Figure 1(b) shows the linear magnetostriction $\Delta L_{i}/L_{i}$
$(i=a,b,{\rm and}~{}c)$ and the volume magnetostriction $\Delta V/V$ at
$T_{\rm ini}$ = 1.4 K for $H~{}||~{}b$. In agreement with the previous study
[39], with increasing magnetic fields $\Delta L_{c}/L_{c}$ increases, while
$\Delta L_{a}/L_{a}$ and $\Delta L_{b}/L_{b}$ decrease at low fields.
Discontinuous changes in the $\Delta L_{i}/L_{i}$s with the $H$-hysteresis
were observed for all the axes, accompanied by the first-order metamagnetic
transition at $H_{\rm m}$ [see also Fig. 1(a)]. Similar to the sign of the
$dL_{i}/dH$ at low fields, an increase in $L_{c}$ and decreases in $L_{a}$ and
$L_{b}$ across $H_{\rm m}$ for the field-up sweep were observed. The largest
discontinuity was observed for the $L_{b}$ as that in the low field region.

Figure 2: (Color online) Pressure dependence of $H_{\rm m}$ of UTe2 for
$H~{}||~{}b$. The data points shown by circles and triangles are taken from
Refs. Knebel2020 and Lin2020, respectively. The dashed line indicates
$\frac{\partial\mu_{0}H_{\rm m}}{\partial P}\approx-11.0$ T/GPa obtained by
the Clausius-Clapeyron’s equation. The solid line is a guide for the eyes.

First, we estimated the volume magnetostriction and compared it with the
reported pressure measurements [31, 32]. Discontinuous change in the volume
across $H_{\rm m}$ at $T_{\rm ini}$ = 1.4 K is evaluated by $\Delta
V/V=\sum_{i=a,b,c}(\Delta L_{i}/L_{i})\approx-5.9\times 10^{-4}$. Adopting
this value and $\Delta M\approx 0.5~{}\mu_{\rm B}$/U [6, 9] to the Clausius-
Clapeyron’s equation, $\frac{\partial\mu_{0}H_{\rm m}}{\partial
p}=\frac{\Delta V}{\Delta M}$, $\frac{\partial\mu_{0}H_{\rm m}}{\partial
p}\approx-11.0~{}{\rm T/GPa}$ is obtained as the initial slope of $H_{\rm m}$
as a function of pressure. Here, we used the reported molar volume $V_{\rm
mol}\approx 5.22\times 10^{-5}~{}{\rm m^{3}/mol}$ at 2.7 K [17]. The direct
comparison of this value to the $p$ dependence of $H_{\rm m}$ is shown in Fig.
2. The $H_{\rm m}(p)$ determined thermodynamically and experimentally agrees
very well [40]. More details on the volume contraction at $H_{\rm m}$ will be
discussed later.

Figure 3: (Color online) (a) $M^{2}$ dependence of $\Delta V/V$ of UTe2 at
$T_{\rm ini}$ = 1.4 K for $H~{}||~{}b$. The red squares (left and bottom axes)
and blue circles (right and top axes) correspond to the data in field regions
below and above $H_{\rm m}$, respectively. There is a kink structure at
$\mu_{0}H=43$ T, indicated by an arrow. The dashed lines indicate a relation,
$\Delta V/V\propto M^{2}$. (b) Field dependence of $\Delta V/V$ at $T_{\rm
ini}$ = 1.4 and 4.2 K. The inset shows $\Delta V/V$ as a function of the
square of the fields at $T_{\rm ini}$ = 1.4 K.

Using the magnetostriction and magnetization results, we can also derive the
pressure dependence of the magnetic susceptibility by employing Maxwell’s
relation, $\left(\frac{\partial M}{\partial
p}\right)_{H,T}=-\left(\frac{\partial V}{\partial\mu_{0}H}\right)_{p,T}$. In a
Pauli paramagnetic system, the magnetic susceptibility $\chi\equiv
dM/d\mu_{0}H$ at low temperature is nearly constant. Thus, the pressure
derivative of the susceptibility is obtained as

$\frac{d\chi}{dp}=-2C_{v}V_{\rm mol}\chi^{2},$ (1)

where $C_{v}$ is the magnetovolume coupling constant, which is obtained by
$\frac{\Delta V}{V}=C_{v}M^{2}$ [41]. From the plot $\Delta V/V$ versus
$M^{2}$ below $H_{\rm m}$ at $T_{\rm ini}$ = 1.4 K shown as red squares in
Fig. 3(a), $C_{v}\approx-3.6\times 10^{-4}~{}(\mu_{\rm B}{\rm/U})^{-2}$ was
evaluated. Using the $C_{v}$ and $\chi$ for $H~{}||~{}b$ [6] to Eq. (1),
$\frac{d\chi}{dp}\approx 1.0\times 10^{-3}$ $\mu_{\rm B}$/(T$\cdot$GPa) is
obtained. This estimation agrees well with the reported pressure dependence of
the magnetic susceptibility [1, 42]. The proportionality of $\Delta V/V\propto
M^{2}$ holds even above $H_{\rm m}$, while a slope changes approximately at 43
T. Although the origin of this anomaly is not clear at present, it may reflect
the change in the 5$f$ characters.

As shown in Fig. 1(b), a strongly anisotropic magnetostriction (AMS) was
observed. Let us see the anisotropy $\Delta L_{i}^{\rm AMS}/L_{i}$ defined by
$\frac{\Delta L_{i}^{\rm AMS}}{L_{i}}\equiv\frac{\Delta
L_{i}}{L_{i}}-\frac{1}{3}\frac{\Delta V}{V}$[43]. The second term is an
averaged lattice change expected for the isotropic case. Figure 1(c) depicts
the AMS along each axis for $H~{}||~{}b$ at $T_{\rm ini}$ = 1.4 K. $\Delta
L^{\rm AMS}_{a}/L_{a}$ shows a weak field dependence, whereas $\Delta L^{\rm
AMS}_{b}/L_{b}$ and $\Delta L^{\rm AMS}_{c}/L_{c}$ show significantly
anisotropic behaviors. At $H_{\rm m}$, a discontinuous rise in $\Delta L^{\rm
AMS}_{c}/L_{c}$ and discontinuous drop in $\Delta L^{\rm AMS}_{b}/L_{b}$ are
more prominent. Interestingly, the absolute values of the AMS for the $b$ and
$c$ axes, including the discontinuous changes at $H_{\rm m}$, are almost
identical. These facts indicate the lattice instability in the $bc$ plane,
which may be responsible for the field-induced SC phases around $H_{\rm m}$
with rotating the field direction between $b$ and $c$ axis [7, 8]. The lattice
instability in the $bc$ plane may also affect the antiferromagnetic
fluctuations with $k_{1}$ = (0, 0.057, 0) [20, 21]. Although the Fermi surface
has not been clarified experimentally [44, 45], the two-dimensionality is
discussed from band calculations [46, 22].

To characterize the electronic state, Grüneisen parameters of a characteristic
energy $X$, $\Gamma_{X}\equiv-\frac{\partial\log{X}}{\partial\log{V}}$, are
useful [47]. Using the volume magnetostriction results below $H_{\rm m}$, the
Grüneisen parameters can be divided into magnetic ($\Gamma_{H}$) and thermal
electronic ($\Gamma_{T}$) contributions [48, 49],

$S_{v}=\frac{1}{2BV_{\rm
mol}}\left[(2\Gamma_{H}-\Gamma_{T})\chi+\Gamma_{T}T\chi^{\prime}\right],$ (2)

where $S_{v}$ is a coefficient of the $H^{2}$-term of $\Delta V/V$, $B=57$ GPa
the bulk modulus [50], and $\chi^{\prime}=\partial\chi/\partial T$. For
$H~{}||~{}b$, more generally in Pauli paramagnets at low temperatures, there
is less temperature dependence of $\chi$. Thus, neglecting the second term of
Eq. (2), we obtain $2\Gamma_{H}-\Gamma_{T}\approx-5$ from our magnetostriction
and magnetization results below $H_{\rm m}$ at $T_{\rm ini}$ = 1.4 K [see the
inset of Fig. 3(b)]: $S_{v}=-5.6\times 10^{-8}~{}{\rm T}^{-2}$ and
$\chi\approx 0.0122~{}{\mu_{\rm B}}/{\rm T}$. Using the reported
$\Gamma_{T}\approx-30$ [39, 51, 1], $\Gamma_{H}\approx-17.5$ is derived.
$\Gamma_{H}$ is nearly half of $\Gamma_{T}$.

Figure 4: (Color online) Comparison between absolute values of the thermal
($\Gamma_{T}$) and magnetic ($\Gamma_{H}$) Grüneisen parameters in various
compounds. The broken line indicates $|\Gamma_{T}|=|\Gamma_{H}|$. As plotted
in different symbols, these intermetallic compounds may be classified into
intermediate valence systems ($|\Gamma_{T}|<|\Gamma_{H}|$, triangles), heavy
fermion systems ($|\Gamma_{T}|\approx|\Gamma_{H}|$, squares), and the others
($|\Gamma_{T}|>|\Gamma_{H}|$, circles). The data except for this study (UTe2)
are obtained from Refs. Kaiser1988, deVisser1992, Dijk1995, Lacerda1989.

Figure 4 compares the relation of the absolute values of $\Gamma_{T}$ and
$\Gamma_{H}$ in various compounds [52, 53, 54, 55]. For many heavy fermion
systems, $\Gamma_{T}=\Gamma_{H}$ is held [47, 52]. Thus, only a single
characteristic energy scale dominates the electronic and the magnetic
properties in many heavy fermion systems. In other words, a dominant
interaction governs the electronic properties in many heavy fermion systems.
Deviations from the correspondence may be seen in the intermediate valence
systems, such as CeSn3, CePd3, and YbCuAl, whose $|\Gamma_{H}|$ is larger than
$|\Gamma_{T}|$. Another example of $|\Gamma_{T}|>|\Gamma_{H}|$ was reported
for unconventional superconductor UBe13[53], where it is well known that the
band structure and non-Fermi liquid state are strongly-field dependent [56,
57]. A sole magnetic property may not link the electronic property. In UTe2,
$\Gamma_{T}$ is nearly double of $\Gamma_{H}$, indicating different channels
of interactions.

From early experimental investigations, the link between $T_{\chi}^{\rm
max}\approx 35$ K and $\mu_{0}H_{\rm m}\approx 35$ T was discussed [5, 6, 31,
9, 1], as in other paramagnetic systems [59]. Recently, a different
characteristic energy scale $T^{\ast}\approx 12$ K was also proposed [60, 61,
62, 51, 63]. A broad peak-like anomaly appeared at $T^{\ast}$ in some
thermodynamic quantities[51]. NMR studies for $H~{}||~{}a$ have revealed some
anomalies at both of $T_{\chi}^{\rm max}$ and $T^{\ast}$; starting a
broadening of NMR spectra below $T_{\chi}^{\rm max}$ and a peak in the nuclear
spin-spin relaxation rate at $T^{\ast}$ [63]. Thus, $T^{\ast}$ may correspond
to the short-range ordering temperature. Interestingly, $\Gamma_{H}$ is almost
identical to $\Gamma_{T}$ at $T^{\ast}$ [51]. This coincidence indicates that
the low field short-range order is eliminated at $H_{\rm m}$.

To further address the influence of the metamagnetic transition on the
electronic properties, we discuss $\Gamma_{T}$ above $H_{\rm m}$. We estimate
the volume thermal expansion using magnetostriction results at $T_{\rm ini}$ =
1.4 and 4.2 K, as shown in Fig. 3(b). Assuming the isothermal process above
$H_{\rm m}$, the volume thermal expansion coefficient $\alpha_{V}\equiv
V^{-1}dV/dT\sim[\Delta V(4.2~{}{\rm K})/V-\Delta V(1.4~{}{\rm
K})/V]/(4.2-1.4)\approx$ $7\times 10^{-6}$ K-1 was estimated at 40 T, slightly
higher than $H_{\rm m}$. Using the field dependence of the electronic
specific-heat coefficient obtained by our previous measurements [24, 9],
$\Gamma_{T}\approx$ 110 at 40 T was obtained[64]. The derived large positive
value of $\Gamma_{T}$ above $H_{\rm m}$ is a mark that the Sommerfeld
coefficient $\gamma\propto C/T$ has a strong maximum at $H_{\rm m}$ [6, 24,
9].

The volume contraction in UTe2 at $H_{\rm m}$ for $H~{}||~{}b$ is strikingly
different from the typical example of heavy-fermion pseudo-metamagnetic
transitions, such as CeRu2Si2 [55] and UPt3 [65], which show an increment in
volume at $H_{\rm m}$. Note that the decrease of $\Delta L_{b}/L_{b}$ suggests
that the uniaxial pressure for the $b$ axis will reduce $H_{\rm m}$. This is
indeed observed at the spin reorientation field $H_{R}$ in URhGe [66, 67]. The
volume contraction of UTe2 at $H_{\rm m}$ may correspond to the loose of the
itinerant character and the resultant valence change.

We discuss a possible link of the metamagnetic transition to the valence
instability. As shown in Fig. 2, $H_{\rm m}$ is suppressed at $p_{c}$.
Increasing pressure from the intermediate valence state (U3+ and U4+) at
$p=0$, the uranium valence increases $\sim$0.1 at $p_{c}$, leading the system
to a more tetravalent 5$f^{2}$ configuration [26, 27]. Thus, a volume
contraction accompanied by the valence change at $p_{c}$ is expected. Such a
contraction is indeed observed across $H_{\rm m}$ at ambient pressure. A
significant change is expected to occur in the 5$f$ uranium magnetism as
functions of pressure and magnetic field.

The observed metamagnetic behaviors can be naturally understood on the basis
of the so-called itinerant-localized duality model for the $5f^{3}$-based
heavy fermion system, in which the $5f^{2}$-electron state with the less
hybridization between conduction and $f$ electrons behaves as $localized$ with
the $5f^{2}$ crystalline electric field level scheme, while the
$5f^{1}$-electron state with the larger hybridization behaves as $itinerant$.
Indeed, in the case of UPd2Al3 (a $5f^{3}$-based heavy fermion system), this
physical picture was shown to work quite well [68]. The crucial point is that
the $itinerant$ $5f^{1}$ component behaves similarly to quasiparticles of
$4f^{1}$-based heavy fermion systems, that is, Ce-based system. UTe2, in the
nearly $5f^{3}$ configuration at an ambient condition, would also exhibit the
duality of the $5f$ electrons, as in UPd2Al3. Particularly, the component of
electrons in the $5f^{1}$-electron state can exhibit a discontinuous
metamagnetic transition associated with a first-order valence transition from
Kondo to mixed valence (MV) state as discussed in Refs. Watanabe2008 and
Watanabe2009 (see also the scheme in Fig. 2). In this case, in the region
below $H_{\rm m}$, 5$f^{1}$ components are in the so-called Kondo state (the
$f$ electron number $n_{f}\simeq$ 1), but in the region above $H_{\rm m}$,
5$f^{1}$ components are in the MV (or valence fluctuating) state. This type of
behavior seems to reproduce the magnetic field and pressure dependence of the
valence of 5$f$ electrons observed in UTe2 across the metamagnetic transition,
as shown in Fig 2. Thus, the electronic correlations decrease drastically with
feedback on the band structure at $H_{\rm m}$. This picture seems to be
supported by the large jump of the magnetoresistance [5, 8], the reduction of
carrier density [71], and the drop of $\gamma$ [24, 9]. The corresponding
changes in Fermi surface and magnetic interactions are ingredients of
metamagnetism. Although the link of valence transition across $H_{\rm m}$ at
$p=0$ and across $p_{c}$ at $H=0$ is not much evident, the same trend of
magnetization-easy-axis switch from the $a$ to the $b$ axes across $H_{\rm m}$
and $p_{c}$ were observed. [6, 9, 1]. Qualitatively it seems a promising road;
quantitatively theoretical progress is now necessary.

In summary, we demonstrated the magnetostriction measurements of UTe2 for
$H~{}||~{}b$. A clear drop in the volume magnetostriction accompanied by the
first-order metamagnetic transition was observed, which satisfactorily agrees
with the pressure dependence of $H_{\rm m}$. Significant anisotropic linear
magnetostriction was also revealed; that for the $a$ axis is almost identical
to the averaged volume magnetostriction, while those for the $b$ and $c$ axes
show larger magnetostriction with a similar magnitude but with opposite sign.
This anisotropy indicates lattice instability within the $bc$ plane, which may
trigger the field-reinforced/reentrant superconductivity. Through the
Grüneisen parameter analyses, we discussed the relation between characteristic
energy scales in temperatures and fields. A possible link of the metamagnetic
transition to a valence instability is considerable on the basis of the dual
nature of the itinerant-localized 5$f$ electrons. This suggests the low- and
high-field electronic and magnetic properties are crucially governed by
different interactions. The magnetic field dependence of the duality deserved
to be clarified.

We thank G. Knebel and W. Knafo for fruitful discussions. This work was
supported by KAKENHI (JP15H05884, JP15H05882, JP15K21732, JP16H04006,
JP15H05745, JP17K05555, JP19H00646, JP20K03854, JP20K20889, JP20J10988,
JP20H00130, and JP20KK0061), ICC-IMR, and ERC starting grant
(NewHeavyFermion).

## References

  * [1] S. Ran ,C. Eckberg,Q.-P. Ding, Y. Furukawa, T. Metz, S.R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Science 365, 684 (2019).
  * [2] D. Aoki, A. Nakamura, F. Honda, D.X. Li,Y. Homma, Y. Shimizu, Y. J. Sato, G. Knebel, J.-P. Brison, A. Pourret, D. Braithwaite, G. Lapertot, Q. Niu, M. Vališka, H. Harima, and J. Flouquet, J. Phys. Soc. Jpn. 88, 043702 (2019).
  * [3] P. F. S. Rosa, A. Weiland, S. S. Fender, B. L. Scott, F. Ronning, J. D. Thompson, E. D. Bauer, and S. M. Thomas, arXiv:2110.06200.
  * [4] D. Aoki, J.-P. Brison, J. Flouquet, K. Ishida, G. Knebel, Y. Tokunaga, Y. Yanase, J. Phys.: Condens. Matter 34, 243002 (2022).
  * [5] W. Knafo, M. Vališka, D. Braithwaite, G. Lapertot, G. Knebel, A. Pourret, J.-P. Brison, J. Flouquet, and D. Aoki, J. Phys. Soc. Jpn. 88, 063705 (2019).
  * [6] A. Miyake, Y. Shimizu, Y. J. Sato, D. Li, A. Nakamura, Y. Homma, F. Honda, J. Flouquet, M. Tokunaga, and D. Aoki, J. Phys. Soc. Jpn. 88, 063706 (2019).
  * [7] S. Ran, I.-L. Liu, Y. S. Eo, D. J. Campbell, P. Neves, W. T. Fuhrman, S. R. Saha, C. Eckberg, H. Kim, J. Paglione, D. Graf, J. Singleton, and N. P. Butch, Nat. Phys. 15, 1250 (2019).
  * [8] W. Knafo, M. Nardone, M. Vališka, A. Zitouni, G. Lapertot, D. Aoki, G. Knebel, and D. Braithwaite, Commun. Phys. 4, 40, (2021).
  * [9] A. Miyake, Y. Shimizu, Y. J. Sato, D. Li, A. Nakamura, Y. Homma, F. Honda, J. Flouquet, M. Tokunaga, and D. Aoki, J. Phys. Soc. Jpn. 90,103702 (2021).
  * [10] G. Knebel, W. Knafo, A. Pourret, Q. Niu, M. Vališka, D. Braithwaite, G. Lapertot, M. Nardone, A. Zitouni, S. Mishra, I. Sheikin, G. Seyfarth, J.-P. Brison, D. Aoki, and J. Flouquet, J. Phys. Soc. Jpn. 88, 063707 (2019).
  * [11] F. Lévy, I. Sheikin, B. Grenier, and A. D. Huxley: Science 309, 1343 (2005).
  * [12] F. Lévy, I. Sheikin, and A. D. Huxley: Nat. Phys. 3, 460 (2007).
  * [13] D. Aoki, T. D. Matsuda, Y. Taufour, E. Hassinger, G. Knebel, and J. Flouquet, J. Phys. Soc. Jpn., 78, 113709 (2009).
  * [14] D. Aoki and J. Flouquet: J. Phys. Soc. Jpn. 81, 011003 (2012).
  * [15] D. Aoki, K. Ishida, and J. Flouquet, J. Phys. Soc. Jpn. 88, 022001 (2019).
  * [16] S. Sundar, S. Gheidi, K. Akintola, A. M. Côté, S.R. Dunsiger, S. Ran, N. P. Butch, S. R. Saha, J. Paglione, and J. E. Sonier, Phys. Rev. B 100, 140502(R) (2019).
  * [17] V. Hutanu, H. Deng, S. Ran, W. T. Fuhrman, H. Thoma, and N. Butch, Acta Cryst. B 76, 137 (2020).
  * [18] C. Paulsen, G. Knebel, G. Lapertot, D. Braithwaite, A. Pourret, D. Aoki, F. Hardy, J. Flouquet, and J.-P. Brison, Phys. Rev. B 103, L180501 (2021).
  * [19] Y. Tokunaga, H. Sakai, S. Kambe, T. Hattori, N. Higa, G. Nakamine, S. Kitagawa, K. Ishida, A. Nakamura, Y. Shimizu, Y. Homma, D. Li, F. Honda, and D. Aoki, J. Phys. Soc. Jpn. 88, 073701 (2019).
  * [20] C. Duan, K. Sasmal, M. B. Maple, A. Podlesnyak, J.-X. Zhu, Q. Si, and P. Dai, Phys. Rev. Lett. 25, 237003 (2020).
  * [21] W. Knafo, G. Knebel, P. Steffens, K. Kaneko, A. Rosuel, J.-P. Brison, J. Flouquet, D. Aoki, G. Lapertot, and S. Raymond, Phy. Rev. B 104, L100409 (2021).
  * [22] Y. Xu, Y. Sheng, and Y.-F. Yang, Phys. Rev. Lett. 123, 217002 (2019).
  * [23] J. Ishizuka and Y. Yanase, Phys. Rev. B 103, 094504 (2021).
  * [24] S. Imajo, Y. Kohama, A. Miyake, C. Dong, J. Flouquet, K. Kindo, and D. Aoki, J. Phys. Soc. Jpn. 88, 083705 (2019).
  * [25] K. Miyake, J. Phys. Soc. Jpn. 90, 024701 (2021).
  * [26] S. M. Thomas, F. B. Santos, M. H. Christensen, T. Asaba, F. Ronning, J. D. Thompson, E. D. Bauer, R. M. Fernandes, G. Fabbris and P. F. S. Rosa, Sci. Adv. 6, eabc8709 (2020).
  * [27] S.-i. Fujimori, I. Kawasaki, Y. Takeda, H. Yamagami, A. Nakamura, Y. Homma, and D. Aoki, J. Phys. Soc. Jpn. 90, 015002 (2021).
  * [28] D. Braithwaite, M. Vališka, G. Knebel, G. Lapertot, J.-P. Brison, A. Pourret, M. E. Zhitomirsky, J. Flouquet, F. Honda, and D. Aoki, Commun. Phys. 2, 147 (2019).
  * [29] S. Ran, H. Kim, I-Lin Liu, S. R. Saha, I. Hayes, T. Metz, Y. S. Eo, J. Paglione, and N. P. Butch, Phys. Rev. B 101, 140503(R) (2020).
  * [30] D. Aoki, F. Honda, G. Knebel, D. Braithwaite, A. Nakamura, D. Li, Y. Homma, Y. Shimizu, Y. J. Sato, J.-P. Brison, and J. Flouquet, J. Phys. Soc. Jpn. 89, 053705 (2020).
  * [31] G. Knebel, M. Kimata, M. Vališka, F. Honda, D. Li, D. Braithwaite, G. Lapertot, W. Knafo, A. Pourret, Y. J. Sato, Y. Shimizu, T. Kihara, J.-P. Brison, J. Flouquet, and D. Aoki, J. Phys. Soc. Jpn. 89, 053707 (2020).
  * [32] W.-C. Lin, D. J. Campbell, S. Ran, I.-L. Liu, H. Kim, A. H. Nevidomskyy, D. Graf, N. P. Butch, and J. Paglione, npj Quantum Mater. 5, 68 (2020).
  * [33] D. Aoki, M. Kimata, Y. J. Sato, G. Knebel, F. Honda, A. Nakamura, D. Li, Y. Homma, Y. Shimizu, W. Knafo, D. Braithwaite, M. Vališka, A. Pourret, J.-P. Brison, and J. Flouquet, J. Phys. Soc. Jpn. 90, 074705 (2021).
  * [34] D. Li, A. Nakamura, F. Honda, Y. J. Sato, Y. Homma, Y. Shimizu, J. Ishizuka, Y. Yanase, G. Knebel, J. Flouquet, and D. Aoki, J. Phys. Soc. Jpn. 90, 073703 (2021).
  * [35] The idiom “linear magnetostriction” used in this manuscript does not mean the magnetostriction in proportion to magnetic fields due to piezomagnetism. We use it to distinguish it from “volume magnetostriction”.
  * [36] A. Ikeda, Y. H. Matsuda, and H. Tsuda, Rev. Sci. Instrum. 89, 096103 (2018).
  * [37] M. Gen, A. Miyake, H. Yagiuchi, Y. Watanabe, A. Ikeda, Y. H. Matsuda, M. Tokunaga, T. Arima, Y. Tokunaga, arXiv:2203.16397.
  * [38] Since the MCE accompanied by the metamagnetic transition is not negligible, even the sample is immersed in the superfluid 4He. The sample temperature above $H_{\rm m}$ changes. [9] Thus, $T_{\rm ini}$ in this manuscript is the temperatures at the beginning of the field generation.
  * [39] S. M. Thomas, C. Stevens, F. B. Santos, S. S. Fender, E. D. Bauer, F. Ronning, J. D. Thompson, A. Huxley, and P. F. S. Rosa, Phys. Rev. B 104, 224501 (2021).
  * [40] The small discrepancy of $H_{\rm m}(p)$ determined thermodynamically and by direct pressure measurements[31, 32] may arise from the uncertainty of the absolute value of $\Delta V/V$ and lack of the $H_{\rm m}(p)$ data near ambient pressure. Coupling between the fiber gauge and sample through glue may cause an insufficient transmission of the sample length change to the gauge. At least in this study, however, we confirmed the reproducibility of $\Delta L_{b}/L_{b}$ with different experimental setups. The sample heating owing to the MCE may also reduce the change in magnetostriction[9]. Considering these uncertainties, the $H_{\rm m}(p)$ determined thermodynamically and experimentally agrees satisfactory. A key lacking measurement is to explore the $H_{\rm m}(p)$ in the low pressure range ($p\leqq 0.1$ GPa).
  * [41] N. V. Mushnikov and T. Goto, Phys. Rev. B 70, 054411 (2004).
  * [42] In the cgs emu units, $\frac{d\chi}{dp}\approx 5.6\times 10^{-4}$ emu/($\rm{mol}\cdot$Oe$\cdot$GPa). The online supplemental material compares $\chi(p)$ determined thermodynamically and experimentally.
  * [43] J. Zieglowski, H. U. Häfner, and D. Wohlleben, Phys. Rev. Lett. 56, 193 (1986).
  * [44] S.-i Fujimori, I. Kawasaki, Y. Takeda, H. Yamagami, A. Nakamura, Y. Homma, and D. Aoki, J. Phys. Soc. Jpn. 88, 103701 (2019).
  * [45] L. Miao, S. Liu, Y. Xu, E. C. Kotta, C. J. Kang, S. Ran, J. Paglione, G. Kotliar, N. P. Butch, J. D. Denlinger, and L. A. Wray, Phys. Rev. Lett. 124, 076401 (2020).
  * [46] J. Ishizuka, S. Sumita, A. Daido, and Y. Yanase, Phys. Rev. Lett. 123, 217001 (2019).
  * [47] J. Flouquet, Prog. Low Temp. Phys. 15, 139 (2005).
  * [48] P. Thalmeier and P. Fulde, Europhys. Lett. 1, 367 (1986).
  * [49] $\Gamma_{H}$ used here is the relative volume dependence of the characteristic magnetic fields. The “magnetic Grüneisen parameter” is often used for a different property, namely the relative field dependence of the characteristic temperature. In such a case, $\Gamma_{H}$ equals the adiabatic magnetocaloric effect.
  * [50] F. Honda, T. Koizumi, N. Kawamura, S. Kobayashi, S. I. Kawaguchi, N. Ishimatsu, Y. J. Sato, Y. Homma, and D. Aoki, private communication.
  * [51] K. Willa, F. Hardy, D. Aoki, D.X. Li, P. Wiecki, G. Lapertot, and C. Meingast, Phys. Rev. B 104, 205107 (2021).
  * [52] A. B. Kaiser and P. Flude, Phys. Rev. B 37, 5357 (1988).
  * [53] A. de Visser, N. H. van Dijk, K. Bakker, J. J. M. Franse, A. Lacerda, J. Flouquet, Z. Fisk, and J. L. Smith, Phys. Rev. B 45, 2962 (1992).
  * [54] N. H. van Dijk, A. de Visser, J. J. M. Franse, and A. A. Menovsky, Phys. Rev. B 51, 12665 (1995).
  * [55] A. Lacerda, A. de Visser, P. Haen, P. Lejay, and J. Flouquet, Phys. Rev. B 40, 8759 (1989).
  * [56] Y. Shimizu, A. Pourret, G. Knebel, A. Palacio-Morales, and D. Aoki, Phys. Rev. B 92, 241101(R) (2015).
  * [57] It is also noted that metamagnetic transition in UBe13 is not observed up to 60 T in Ref. Detwiler2000.
  * [58] J. A. Detwiler, G. M. Schmiedeshoff, N. Harrison, A. H. Lacerda, J. C. Cooley, and J. L. Smith, Phys. Rev. B 61, 402 (2000).
  * [59] D. Aoki, W. Knafo, and I. Sheikin, C. R. Phys. 14, 53 (2013).
  * [60] Q. Niu, G. Knebel, D. Braithwaite, D. Aoki, G. Lapertot, G. Seyfarth, J.-P. Brison, J. Flouquet, and A. Pourret, Phys. Rev. Lett. 124, 086601 (2020).
  * [61] L. P. Cairns, C. R. Stevens, C. D. O’Neill, and A. Huxley, J. Phys.: Condens. Matter 32, 415602 (2020).
  * [62] Y. S. Eo, S. R. Saha, H. Kim, S. Ran, J. A. Horn, H. Hodovanets, J. Collini, W. T. Fuhrman, A. H. Nevidomskyy, N. P. Butch, M. S. Fuhrer, and J. Paglione, arXiv:2101.03102.
  * [63] Y. Tokunaga, H. Sakai, S. Kambe, Y. Haga, Y. Tokiwa, P. Opletal, H. Fujibayashi, K. Kinjo, S. Kitagawa, K. Ishida, A. Nakamura, Y. Shimizu, Y. Homma, D. Li, F. Honda, and D. Aoki, J. Phys. Soc. Jpn. 91, 023707 (2022).
  * [64] Here, we follow the conventional definition of the Grüneisen parameter, $\Gamma_{T}=\frac{V_{\rm mol}\alpha_{V}}{\kappa C}$, where $\kappa=1/B$ and $C$ are compressibility and specific heat, respectively.
  * [65] A. de Visser, L. Puech, W. Joss, A. A. Menovsky, and J. J. M. Franse, Jpn. J. Appl. Phys. 26, 513 (1987).
  * [66] D. Braithwaite, D. Aoki, J.-P. Brison, J. Flouquet, G. Knebel, A. Nakamura, and A. Pourret, Phys. Rev. Lett. 120, 037001 (2018).
  * [67] S. Nakamura, T. Sakakibara, Y. Shimizu, S. Kittaka, Y. Kono, Y. Haga, J. Pospíšil, and E. Yamamoto, Prog. Nucl. Sci. Tech. 5, 123 (2018).
  * [68] N. K. Sato, N. Aso, K. Miyake, R. Shiina, P. Thalmeier, G. Varelogiannis, C. Geibel, F. Steglich, P. Fulde, and T. Komatsubara, Nature 410, 340 (2001).
  * [69] S. Watanabe, A. Tsuruta, K. Miyake, and J. Flouquet, Phys. Rev. Lett. 100, 236401 (2008).
  * [70] S. Watanabe, A. Tsuruta, K. Miyake, and J. Flouquet, J. Phys. Soc. Jpn. 78, 104706 (2009).
  * [71] Q. Niu, G. Knebel, D. Braithwaite, D. Aoki, G. Lapertot, M. Vališka, G. Seyfarth, W. Knafo, T. Helm, J.-P. Brison, J. Flouquet, and A. Pourret, Phys. Rev. Reserch 2, 033179 (2020).

## Supplemental material

##  1 Comparison of $\chi(p)$ determined thermodynamically and experimentally

We compare the pressure dependence of the magnetic susceptibility $\chi$ for
$H~{}||~{}b$ determined thermodynamically and experimentally. Figure 5 shows
pressure dependence of $\chi_{0}$ at 10000 Oe for $H~{}||~{}b$, which is
estimated by extrapolating to 0 K [1]. From Maxwell’s relation using our
magnetization and magnetostriction results, $\frac{d\chi}{dp}=5.6\times
10^{-4}$ emu/(mol$\cdot$ Oe$\cdot$ GPa) is obtained. As shown by a dashed line
in Fig. 5, our evaluated $\chi_{0}(p)$ agrees nicely with pressure dependence
of $\chi_{0}$.

Figure 5: Pressure dependence of $\chi_{0}$ at 10000 Oe for $H~{}||~{}b$ of
UTe2. The dashed line indicates $d\chi/dp\approx 5.6\times 10^{-4}$
emu/(mol$\cdot$ Oe$\cdot$ GPa). The vertical dashed line indicates a critical
pressure $p_{c}$. The $\chi_{0}(p)$ data are taken from Ref. Li2021.

## References

  * [1] D. Li, A. Nakamura, F. Honda, Y. J. Sato, Y. Homma, Y. Shimizu, J. Ishizuka, Y. Yanase, G. Knebel, J. Flouquet, and D. Aoki, J. Phys. Soc. Jpn. 90, 073703 (2021).


#  Sensing With Random Signals

###### Abstract

Radar systems typically employ well-designed deterministic signals for target
sensing. In contrast to that, integrated sensing and communications (ISAC)
systems have to use random signals to convey useful information, potentially
causing sensing performance degradation. This paper analyzes the sensing
performance via random ISAC signals over a multi-antenna system. Towards this
end, we define a new sensing performance metric, namely, ergodic linear
minimum mean square error (ELMMSE), which characterizes the estimation error
averaged over the randomness of ISAC signals. Then, we investigate a data-
dependent precoding scheme to minimize the ELMMSE, which attains the optimized
sensing performance at the price of high computational complexity. To reduce
the complexity, we present an alternative data-independent precoding scheme
and propose a stochastic gradient projection (SGP) algorithm for ELMMSE
minimization, which can be trained offline by locally generated signal
samples. Finally, we demonstrate the superiority of the proposed methods by
simulations.

Index Terms—  Integrated sensing and communications (ISAC), deterministic-
random tradeoff, Gaussian signals.

##  1 Introduction

Recently, the International Telecommunication Union (ITU) has recognized
integrated sensing and communications (ISAC) as one of the six key usage
scenarios for 6G, marking the continuous progression of ISAC from concept to
practice [1]. By reusing wireless resources, ISAC may deploy ubiquitous
sensing ability over existing networking infrastructures in a low-cost manner
[2, 3, 4]. Among various approaches to implement ISAC, utilizing communication
signals for sensing holds the most promising potential [5, 6, 7]. To convey
communication information, ISAC signals should be random, e.g., capacity-
achieving Gaussian signals [8]. Conversely, radar systems prefer deterministic
signals, e.g., constant-modulus signals, which may offer improved sensing
performance. This contradictory conflict leads to the so-called deterministic-
random tradeoff (DRT) in ISAC systems [9].

Despite the random nature of ISAC signals, most of previous studies focusing
on the ISAC signal design treat the sample covariance matrix of transmit
signals as deterministic [10, 11, 12], under the assumption that the data
frame is sufficiently long. Those approaches overlooked the DRT, making them
akin to infinite-frame-length approximations. However, such approximations may
not be reliable in practical scenarios. For instance, in massive multiple-
input multiple-output (MIMO) systems, the base station may be equipped with
massive antennas, and the signal sample covariance matrix converges to its
statistical counterpart only when the frame length is much larger than the
antenna number. To alleviate the resultant huge computational and signal
processing overheads, one has to perform sensing by leveraging short frames,
where the law of large numbers no longer holds and the data randomness cannot
be eliminated. This motivates our investigation in this paper.

Specifically, we consider sensing with random signals in MIMO ISAC systems,
where the transmitter aims to estimate the target impulse response (TIR)
matrix by using random information-carrying communication signals. To
characterize the sensing performance in this practical scenario, we define a
new performance metric, namely, ergodic linear minimum mean-square-error
(ELMMSE), which is defined as the average estimation error over random signal
realizations. Then, we reveal the performance disparity between deterministic
and random signals for sensing, drawing on Jensen’s inequality. Subsequently,
we present a data-dependent precoding strategy and a stochastic gradient
projection (SGP) algorithm to minimize the ELMMSE. Finally, our simulations
underscore the performance degradation linked to random signals in sensing,
and offer guidance for practical design in scenarios with limited coherent
frame lengths.

##  2 System Modeling

As illustrated in Fig. 1, we consider a MIMO ISAC system with the base station
(BS) deployed with $N_{t}$ transmit antennas and $N_{r}$ receive antennas,
where the target sensing is implemented over a coherent processing interval
consisting of $L$ snapshots in a monostatic setup. The echo signals at the
ISAC receiver are denoted by

$\displaystyle\bm{Y}_{s}$ $\displaystyle=\bm{H}_{s}\bm{X}+\bm{Z}_{s},$ (1)

where $\bm{X}\in\mathbb{C}^{N_{t}\times L}$ is the sensing signal matrix,
$\bm{H}_{s}\in\mathbb{C}^{N_{r}\times N_{t}}$ is the TIR matrix to be
estimated, and $\bm{Z}_{s}\in\mathbb{C}^{N_{r}\times L}$ denotes an additive
noise matrix with each entry following $\mathcal{CN}(0,\sigma_{s}^{2})$. In
what follows, we commence by sensing with deterministic signals.

Fig. 1: The illustration of considered MIMO ISAC system using Gaussian
signals.

###  2.1 Sensing With Deterministic Training Signals

Let $\bm{S}_{D}\in\mathbb{C}^{N_{t}\times L}$ denote a deterministic training
signal satisfying $(1/L)\bm{S}_{D}\bm{S}_{D}^{H}=\bm{I}_{N_{t}}$. The sensing
signal matrix is $\bm{X}=\bm{W}\bm{S}_{D}$, where
$\bm{W}\in\mathbb{C}^{N_{t}\times N_{t}}$ represents the precoding matrix to
be optimized. Then, the celebrated linear minimum mean-square-error (LMMSE)
estimation of $\bm{H}_{s}$ is

$\displaystyle\hat{\bm{H}}_{\mathsf{LMMSE}}=\bm{Y}_{s}\left(\bm{X}^{H}\bm{R}_{H}\bm{X}+\sigma_{s}^{2}N_{r}\bm{I}_{L}\right)^{-1}\bm{X}^{H}\bm{R}_{H},$
(2)

where $\bm{R}_{H}=\mathbb{E}[\bm{H}_{s}^{H}\bm{H}_{s}]$ represents the
statistical correlation matrix of the channel with $\mathbb{E}[\cdot]$
denoting the statistic expectation. The resulting estimation error is [13]

$\displaystyle{J}_{\mathsf{LMMSE}}=\mathrm{tr}\Big{[}\Big{(}\bm{R}_{H}^{-1}+\frac{L}{\sigma_{s}^{2}N_{r}}\bm{W}\bm{W}^{H}\Big{)}^{-1}\Big{]}.$
(3)

Accordingly, the LMMSE-oriented design is to solve the deterministic
optimization problem:

$\displaystyle\min_{\bm{W}}~{}J_{\mathsf{LMMSE}}~{}~{}\mathrm{s.t.}~{}\bm{W}\in\mathcal{A},$
(4)

where $\mathcal{A}=\\{\bm{W}\in\mathbb{C}^{N_{t}\times
N_{t}}|~{}\|\bm{W}\|_{F}^{2}\leq P\\}$ denotes the feasible region and $P$ is
the power budget. Let $\bm{R}_{H}=\bm{Q}\bm{\varLambda}\bm{Q}^{H}$ denote the
eigenvalue decomposition of $\bm{R}_{H}$. The optimal precoding matrix is
known to be the following water-filling solution [13]:

$\displaystyle\bm{W}_{\mathsf{LMMSE}}^{\mathsf{opt}}=\sqrt{\frac{\sigma_{s}^{2}{N_{r}}}{L}}\bm{Q}\Big{[}\Big{(}\mu_{0}\bm{I}_{N_{t}}-\bm{\varLambda}^{-1}\Big{)}^{+}\Big{]}^{\frac{1}{2}},$
(5)

where $(x)^{+}=\max(x,0)$ and $\mu_{0}$ is a constant that is determined based
on $\|\bm{W}_{\mathsf{LMMSE}}^{\mathsf{opt}}\|_{F}^{2}=P$.

###  2.2 Sensing With Random Communication Signals

In contrast to conventional radar systems, ISAC systems have to employ random
communication signals for sensing. In this paper, we consider Gaussian
signaling for ISAC systems, which is known to achieve the capacity of point-
to-point Gaussian channels [8]. Let
$\bm{S}=[\bm{s}(1),\dots,\bm{s}(L)]\in\mathbb{C}^{N_{t}\times L}$ denote the
transmitted ISAC signal, where each column is independent and identically
distributed (i.i.d.) and follows the circularly symmetric complex Gaussian
(CSCG) distribution with zero mean and covariance $\bm{I}_{N_{t}}$, i.e.,
$\bm{s}(l)\sim\mathcal{CN}(\mathbf{0},\bm{I}_{N_{t}})$. It is worth noting
that the following assumption is commonly seen in the existing works [10, 11,
12],

$\displaystyle\frac{1}{L}\bm{S}\bm{S}^{H}\approx\mathbb{E}[\bm{s}(l)\bm{s}^{H}(l)]=\bm{I}_{N_{t}},l=1,\dots,L,$
(6)

where $L\gg N_{t}$ is assumed such that the approximation holds. However, this
assumption cannot always be satisfied, especially when the BS does not have
the ability to accumulate sufficiently long frames. To this end, one needs to
take the randomness of the Gaussian signals into account when designing the
precoding matrix. More severely, the sensing performance may be randomly
varying due to the random signaling. As a consequence, it may not be
appropriate to use the conventional estimation metric relying on the
instantaneous realization of signals. Therefore, we define a new ELMMSE to
measure the average sensing performance over random signal realizations, which
is expressed as

$\displaystyle{J}_{\mathsf{ELMMSE}}=\mathbb{E}_{\bm{S}}\left\\{\mathrm{tr}\Big{[}(\bm{R}_{H}^{-1}+\frac{1}{\sigma_{s}^{2}N_{r}}\bm{W}\bm{S}\bm{S}^{H}\bm{W}^{H})^{-1}\Big{]}\right\\}.$
(7)

In this paper, we aim to minimize the average sensing performance under
Gaussian signaling, which is modeled by 111 In general, $\bm{W}$ can be
designed depending on random signal realization $\bm{S}$. This thus leads to
two different designs of $\bm{W}$, which remains unchanged and changes
adaptively based on the instantaneous value of $\bm{S}$, referred to as data-
independent and data-dependent designs in our paper, respectively.

$\displaystyle\min_{\bm{W}}~{}~{}f(\bm{W})\triangleq
J_{\mathsf{ELMMSE}}~{}~{}\mathrm{s.t.}~{}\bm{W}\in\mathcal{A}.$ (8)

Proposition 1: Gaussian signals lead to ELMMSE that is no lower than the LMMSE
based on the deterministic signals, i.e., LMMSE serves as a lower bound of
ELMMSE.

Proof: By applying the Jensen’s inequality to ELMMSE, we have the following
relationship

$\displaystyle{J}_{\mathsf{ELMMSE}}$
$\displaystyle\overset{(a)}{\geq}\mathrm{tr}\Big{[}\Big{(}\bm{R}_{H}^{-1}+\frac{1}{\sigma_{s}^{2}N_{r}}\mathbb{E}_{\bm{S}}\Big{\\{}\bm{W}\bm{S}\bm{S}^{H}\bm{W}^{H}\Big{\\}}\Big{)}^{-1}\Big{]}$
$\displaystyle\overset{(b)}{=}\mathrm{tr}\Big{[}\Big{(}\bm{R}_{H}^{-1}+\frac{L}{\sigma_{s}^{2}N_{r}}\bm{W}\bm{W}^{H}\Big{)}^{-1}\Big{]}={J}_{\mathsf{LMMSE}},$
(9)

where $(a)$ holds due to the convexity of ${J}_{\mathsf{ELMMSE}}$ with respect
to $\bm{S}\bm{S}^{H}$, and $(b)$ holds following the fact that
$\mathbb{E}_{\bm{S}}[\bm{S}\bm{S}^{H}]=L\bm{I}_{N_{t}}$ [14]. This completes
the proof. $\square$

Remark 1: We notice that the Jensen’s inequality in (2.2) tends to be tight
with $L$ going to infinity. In light of Proposition 1, we show the asymptotic
performance of random signals in Fig. 2, precoded by the water-filling
precoder given in (5), indicating that the assumption in (6) is not always
reliable. It also justifies the use of ELMMSE instead of LMMSE as the
performance metric or design objective under random signaling, since
minimizing the LMMSE only minimizes the lower bound of the ELMMSE instead of
the exact value. This observation motivates us to present efficient ELMMSE-
oriented precoding solutions for solving problem (8) in the next section.

Fig. 2: The asymptotic performance of Gaussian signals.

##  3 ELMMSE-Oriented precoding

This section presents efficient ELMMSE-oriented precoding solutions for
solving problem (8), namely, data-dependent precoding and SGP-based precoding,
as detailed below.

###  3.1 Data-Dependent Precoding

We first take the stochastic approximation into account by $N$ Gaussian signal
samples, written as

$\displaystyle{f}(\bm{W})\approx\frac{1}{N}\sum_{n=0}^{N-1}\hat{f}(\bm{W};\bm{S}_{n}),$
(10)

where
$\hat{f}(\bm{W};\bm{S}_{n})=\mathrm{tr}[(\bm{R}_{H}^{-1}+\frac{1}{\sigma_{s}^{2}N_{r}}\bm{W}\bm{S}_{n}\bm{S}_{n}^{H}\bm{W}^{H})^{-1}]$
and $\bm{S}_{n}$ denotes the $n$-th Gaussian data realization. Observe that
each $\bm{S}_{n}$ is known to both the transmitter and sensing receiver due to
the monostatic sensing setup. Then the precoding matrix $\bm{W}$ can be seen
as a function of $\bm{S}$ over different data realization indices $n$, written
as $\bm{W}(\bm{S}_{n})$. Therefore, based on (10), problem (8) is approximated
as

$\displaystyle\min_{\\{\bm{W}(\bm{S}_{n})\\}_{n=1}^{N}}~{}\frac{1}{N}\sum_{n=0}^{N-1}\hat{f}(\bm{W}(\bm{S}_{n});\bm{S}_{n})~{}~{}\mathrm{s.t.}~{}\bm{W}(\bm{S}_{n})\in\mathcal{A}.$
(11)

Consequently, the precoding matrix can be optimized successively with given
$\bm{S}_{n}$ over different data realization $n$, which is called data-
dependent precoding. Then this problem is decomposed following $N$ parallel
deterministic sub-problems:

$\displaystyle\min_{\bm{W}_{n}}~{}\hat{f}(\bm{W}_{n};\bm{S}_{n})~{}~{}\mathrm{s.t.}~{}\bm{W}_{n}\in\mathcal{A},$
(12)

where we define $\bm{W}_{n}\triangleq\bm{W}(\bm{S}_{n})$ for notational
convenience. Although problem (12) is non-convex, we propose to solve it using
the successive convex approximation (SCA) algorithm. Let us first approximate
the objective function by its first-order Taylor expansion as

$\displaystyle\hat{f}(\bm{W}_{n};\bm{S}_{n})\approx\hat{f}($
$\displaystyle\bm{W}_{n}^{(t)};\bm{S}_{n})$
$\displaystyle+\langle\nabla{\hat{f}}(\bm{W}_{n}^{(t)};\bm{S}_{n}),(\bm{W}_{n}-\bm{W}_{n}^{(t)})\rangle,$
(13)

where $\langle\bm{A},\bm{B}\rangle$ represents the inner product of $\bm{A}$
and $\bm{B}$, $\bm{W}_{n}^{(t)}\in\mathcal{A}$ denotes a given feasible point,
and $\nabla{\hat{f}}(\cdot)$ represents the gradient of $\hat{f}(\cdot)$,
calculated by

$\displaystyle\nabla\hat{f}(\bm{W}_{n}^{(t)};\bm{S}_{n})=-\frac{1}{\sigma_{s}^{2}N_{r}}\bm{A}_{n}^{-2}\bm{W}_{n}^{(t)}\bm{S}_{n}\bm{S}_{n}^{H},$
(14)

where
$\bm{A}_{n}=\bm{R}_{H}^{-1}+\frac{1}{\sigma_{s}^{2}N_{r}}\bm{W}_{n}^{(t)}\bm{S}_{n}\bm{S}_{n}^{H}(\bm{W}_{n}^{(t)})^{H}$.

At the $t$-th iteration of the SCA algorithm, we solve the following
optimization problem:

$\displaystyle\min_{\bm{W}_{n}}$
$\displaystyle~{}~{}{g}(\bm{W}_{n})\triangleq\langle\nabla{\hat{f}}(\bm{W}_{n}^{(t)};\bm{S}_{n}),(\bm{W}_{n}-\bm{W}_{n}^{(t)})\rangle$
$\displaystyle\mathrm{s.t.}$ $\displaystyle~{}~{}\bm{W}_{n}\in\mathcal{A}.$
(15)

We observe that problem (3.1) is convex and can be solved by off-the-shelf
numerical tools like CVX [15]. We denote the solution to problem (3.1) as
$\bm{W}_{n}^{\prime}\in\mathcal{A}$. We observe that
${g}(\bm{W}_{n}^{\prime})\leq 0$, since ${g}(\bm{W}_{n}^{(t)})=0$, which
indicates that $(\bm{W}_{n}^{\prime}-\bm{W}_{n}^{(t)})$ represents a descent
direction for the $(t+1)$-th iteration. Therefore, we can update the precoding
matrix by moving along this direction with a certain step size $\delta^{(t)}$,
i.e.,

$\displaystyle\bm{W}_{n}^{(t+1)}=\bm{W}_{n}^{(t)}+\delta^{(t)}(\bm{W}_{n}^{\prime}-\bm{W}_{n}^{(t)}),~{}\delta^{(t)}\in[0,1].$
(16)

As detailed in Algorithm 1, by solving a series of parallel sub-problems (12),
we obtain the set of data-dependent precoding matrix tailored for each
$\bm{S}_{n}$, denoted by

$\displaystyle\mathcal{W}=\\{\bm{W}_{1}^{\star},\dots,\bm{W}_{N}^{\star}\\}.$
(17)

Accordingly, the ELMMSE may be calculated as

$\displaystyle{J}_{\mathsf{ELMMSE}}=\frac{1}{N}\sum_{n=0}^{N-1}\hat{f}(\bm{W}^{\star}_{n},\bm{S}_{n}).$
(18)

Remark 2: The data-dependent precoding represents the performance lower bound
of (11). Notice that $\mathcal{W}$ needs to be designed sequentially for
different realizations of transmitted Gaussian signals. This, however, results
in high computational complexity in practical scenarios. This motivates us to
harness the popular SGP algorithm with lower complexity, where a single data-
independent precoding matrix is employed for all data realizations, as
detailed below.

Algorithm 1 SCA Algorithm for Solving (12).

1:$\bm{S}_{n},\bm{R}_{H},P,N,\sigma^{2}_{s},N_{t},N_{r},L,t_{\mathsf{max}},\xi<0$.

2:$\bm{W}^{\star}_{n}=\bm{W}_{n}^{(t)}$.

3:Initialize the precoding matrix $\bm{W}_{n}^{(t)},r=t$.

4:Initialize the Gaussian signal set
$\mathcal{S}=\\{\bm{S}_{1},\ldots\,\bm{S}_{N}\\}$.

5:repeat

6: Solve the convex Problem (3.1) to obatin
$\bm{W}_{n}^{\prime}\in\mathcal{A}$ .

7: Update
$\bm{W}_{n}^{(t+1)}=\bm{W}^{(t)}+\delta^{(t)}(\bm{W}_{n}^{\prime}-\bm{W}^{(t)})$,
where the step $\delta^{t}$ can be set by applying the exact line search.

8: Update $t=t+1$.

9:until ${g}(\bm{W}_{n}^{(t)})\leq\xi$ or $t=t_{\mathsf{max}}$.

###  3.2 SGP-Based Data-Independent Precoding

To seek for a data-independent precoder, we resort to the stochastic gradient
descent (SGD) algorithm, which yields a solution with only one or mini-batch
samples in each iteration, reducing computational complexity. At the $r$-th
iteration, the precoding matrix $\bm{W}$ is updated as

$\displaystyle\bm{W}^{(r+1)}\leftarrow\bm{W}^{(r)}-\eta^{(r)}\hat{\nabla}{f}(\bm{W}^{(r)}),$
(19)

where $\bm{W}^{(r)}$ denotes the precoding matrix in iteration $r$,
$\eta^{(r)}$ represents the step size, and $\hat{\nabla}{f}(\bm{W}^{(r)})$ is
the gradient obtained from a local randomly generated mini-batch Gaussian
sample set $\mathcal{D}^{(r)}$, calculated by

$\displaystyle\hat{\nabla}{f}(\bm{W}^{(r)})=\frac{1}{|\mathcal{D}^{(r)}|}\sum_{\bm{S}\in\mathcal{D}^{(r)}}{\nabla}\hat{f}(\bm{W}^{(r)};\bm{S}).$
(20)

In order to adhere to the transmit power budget, we propose to utilize the SGP
algorithm, where we project $\bm{W}^{(r+1)}$ onto the convex feasible region
$\mathcal{A}$, i.e.,

$\displaystyle\mathsf{Proj}_{\mathcal{A}}(\bm{W})=\left\\{\begin{array}[]{l}\bm{W},~{}\mathrm{if}~{}\bm{W}\in\mathcal{A},\\\
\bm{W}\sqrt{\frac{P}{\|\bm{W}\|_{F}^{2}}},~{}\mathrm{Otherwise}.\end{array}\right.$
(23)

Following the iteration format in (19) and the projection step in (23), we are
now ready to introduce the proposed SGP algorithm to solve (8), as detailed in
Algorithm 2.

Algorithm 2 SGP Algorithm for Solving (8).

1:$\bm{R}_{H},P,N,\sigma^{2}_{s},N_{t},N_{r},r_{\mathrm{max}},\epsilon$.

2:$\bm{W}$.

3:Initialize the precoding matrix $\bm{W}^{(r)},r=1$.

4:Initialize the Gaussian signal set
$\mathcal{S}=\\{\bm{S}_{1},\ldots\,\bm{S}_{N}\\}$.

5:repeat

6: Generate $\mathcal{D}^{(r)}$ and calculate
$\hat{\nabla}{f}(\bm{W}^{(r)},\bm{S}^{(r)})$.

7: Update
$\bm{W}^{(r+1)}=\mathsf{Proj}_{\mathcal{A}}\Big{(}\bm{W}^{(r)}-\eta^{(r)}\hat{\nabla}{f}(\bm{W}^{r})\Big{)}$.

8: Update $r=r+1$.

9:until The increase of the objective value is below $\epsilon$ or
$r=r_{\mathrm{max}}$.

Table 1: Parameters in Simulations. Parameter |  Value |  Parameter |  Value | Parameter |  Value
---|---|---|---|---|---
$N_{t}$ |  $64$ |  $N_{r}$ |  $32$ | $N$ |  $100$
$\sigma_{s}^{2}$ |  $0~{}\mathrm{dBm}$ |  $P$ |  $30~{}\mathrm{dBm}$ | $\xi$ |  $-0.1$
$t_{\mathsf{max}}$ |  $30$ |  $r_{\mathsf{max}}$ |  $2000$ | $\epsilon$ |  $10^{-5}$

##  4 Simulation Results

In this section, we demonstrate the numerical results of the proposed methods.
If not otherwise specified, the transmit SNR is defined as
$LP/\sigma_{s}^{2}$. The simulation parameters are listed in Table 1. For the
SGP method, we set the number of mini-batch samples as
$|\mathcal{D}^{(r)}|=10$. The eigenvalues of $\bm{R}_{H}$ represent the
spatial channel correlation, which follow a uniform distribution on the
interval $[1,10]$ in the simulation.

(a) SCA Algorithm.

(b) SGP Algorithm.

Fig. 3: Convergence examples of proposed algorithms.

(a) $L=32,N_{t}=64$.

(b) $L=64,N_{t}=64$.

Fig. 4: Performance comparison among different schemes versus transmit SNR.

In Fig. 3, we set the frame length $L=32$ and transmit SNR $=30~{}\mathrm{dB}$
to show the convergence of our proposed algorithms in Section 3. In Fig. 3(a),
we illustrate the convergence of our proposed SCA algorithm for one of the
parallel sub-problems (12) within 30 iterations. In Fig. 3(b), we show the
convergence of our proposed SGP algorithm, which can be implemented offline
based on the locally generated Gaussian signal samples, and the step size is
set as $\eta^{(r)}=10/(10+r)$ to guarantee the convergence [16]. It is
observed that SGP converges quickly within tens of iterations.

In Fig. 4, we compare three precoding schemes, namely, classical water-
filling, SGP, and data-dependent precoding. Notably, the SGP algorithm
outperforms the classical water-filling scheme in both cases with $L=32$ and
$L=64$. This is attributed to two key factors. When $L=32$, the rank-deficient
nature of the signal (see Fig. 4(a)) causes an error floor in the water-
filling scheme, rendering it suboptimal even for deterministic signals.
Secondly, the inherent randomness of signals in ISAC scenarios is disregarded
by the classical water-filling scheme. Moreover, the data-dependent precoding
achieves the best performance, leveraging precise knowledge of transmitted
data. Nevertheless, adapting data-dependent precoding for varying data
realizations requires complex redesigns, incurring significant complexity.
Fortunately, a data-independent precoder can be implemented offline through
the SGP, with a favorable performance-complexity tradeoff.

Fig. 5: Performance comparison between deterministic and random signals versus
the transmit SNR.

In Fig. 5, we show the resultant estimation error by transmitting Gaussian
(communication) and deterministic (training) signals, with an increasing SNR.
All the curves of random signals are obtained by SGP. It is observed that with
the increasing frame length, the average sensing performance is asymptotically
close to the LMMSE using deterministic signals since the randomness of the
signals is decreased. However, there is still a performance gap between using
deterministic and random signals, already represented by the Jensen’s
inequality.

##  5 Conclusion

In this paper, we introduced a new performance metric, namely, ELMMSE, to
assess the performance of employing random signals for target sensing. We
first revealed that there is a certain performance gap between deterministic
and random signals used in sensing, a revelation enabled by Jensen’s
inequality. Then, we shed light on the precoding matrix design while taking
the randomness of the transmit signals into account. To minimize the
predefined ELMMSE, we proposed a data-dependent precoding approach and a data-
independent SGP scheme, respectively. Finally, simulation results indicated
the superiority of our proposed design schemes. In ISAC systems, we ascertain
the necessity of accounting for randomness when utilizing Gaussian signals for
sensing, particularly for scenarios with short frame length.

## References

  * [1] I.-R. WP5D, “Draft new Recommendation ITU-R M.[IMT.FRAMEWORK FOR 2030 AND BEYOND],” 2023.
  * [2] Y. Cui, F. Liu, X. Jing, and J. Mu, “Integrating sensing and communications for ubiquitous IoT: Applications, trends, and challenges,” _IEEE Network_ , vol. 35, no. 5, pp. 158–167, 2021.
  * [3] A. Liu, Z. Huang, M. Li, Y. Wan, W. Li, T. X. Han, C. Liu, R. Du, D. K. P. Tan, J. Lu _et al._ , “A survey on fundamental limits of integrated sensing and communication,” _IEEE Commun. Surveys Tuts._ , 2022.
  * [4] S. Lu, F. Liu, Y. Li, K. Zhang, H. Huang, J. Zou, X. Li, Y. Dong, F. Dong, J. Zhu _et al._ , “Integrated sensing and communications: Recent advances and ten open challenges,” _arXiv preprint arXiv:2305.00179_ , 2023\.
  * [5] C. Sturm and W. Wiesbeck, “Waveform design and signal processing aspects for fusion of wireless communications and radar sensing,” _Proceedings of the IEEE_ , vol. 99, no. 7, pp. 1236–1259, 2011.
  * [6] X. Li, V. C. Andrei, U. J. Mönich, and H. Boche, “Optimal linear precoder design for MIMO-OFDM integrated sensing and communications based on Bayesian Cramér-Rao bound,” _arXiv preprint arXiv:2308.12106_ , 2023\.
  * [7] Y. Zhang, S. Aditya, and B. Clerckx, “Input distribution optimization in OFDM dual-function radar-communication systems,” _arXiv preprint arXiv:2305.06635_ , 2023.
  * [8] T. M. Cover, _Elements of Information Theory_. John Wiley & Sons, 1999.
  * [9] Y. Xiong, F. Liu, Y. Cui, W. Yuan, T. X. Han, and G. Caire, “On the fundamental tradeoff of integrated sensing and communications under Gaussian channels,” _IEEE Trans. Inf. Theory_ , vol. 69, no. 9, pp. 5723–5751, 2023.
  * [10] F. Liu, Y.-F. Liu, A. Li, C. Masouros, and Y. C. Eldar, “Cramér-Rao bound optimization for joint radar-communication beamforming,” _IEEE Trans. Signal Process._ , vol. 70, pp. 240–253, 2021.
  * [11] F. Liu, Y.-F. Liu, C. Masouros, A. Li, and Y. C. Eldar, “A joint radar-communication precoding design based on Cramér-Rao bound optimization,” in _2022 IEEE Radar Conference (RadarConf22)_. IEEE, 2022, pp. 1–6.
  * [12] S. Lu, X. Meng, Z. Du, Y. Xiong, and F. Liu, “On the performance gain of integrated sensing and communications: A subspace correlation perspective,” in _2023 IEEE International Conference on Communications (ICC)_. IEEE, 2023.
  * [13] M. Biguesh and A. B. Gershman, “Training-based MIMO channel estimation: A study of estimator tradeoffs and optimal training signals,” _IEEE Trans. Signal Process._ , vol. 54, no. 3, pp. 884–893, 2006.
  * [14] D. Maiwald and D. Kraus, “On moments of complex wishart and complex inverse Wishart distributed matrices,” in _1997 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , vol. 5, 1997, pp. 3817–3820.
  * [15] S. P. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.
  * [16] L. eon Bottou, “Online learning and stochastic approximations,” _Online Learning in Neural Networks_ , vol. 17, no. 9, p. 142, 1998.


# Machine Learning for the LHCb Simulation

Lucio Anderlini

###### Abstract

Most of the computing resources pledged to the LHCb experiment at CERN are
necessary to the production of simulated samples used to predict resolution
functions on the reconstructed quantities and the reconstruction and selection
efficiency. Projecting the Simulation requests to the years following the
upcoming LHCb Upgrade, the relative computing resources would exceed the
pledges by more than a factor of 2. In this contribution, I discuss how
Machine Learning can help to speed up the Detector Simulation for the upcoming
Runs of the LHCb experiment.

The LHCb detector [1, 2] is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, designed for the study of particles
containing $b$ and $c$ quarks. It includes a silicon-strip vertex detector
surrounding the $pp$ interaction region that allows $c$ and $b$ hadrons to be
identified from their characteristically long flight distance; a tracking
system that provides a measurement of momentum, $p$, of charged particles. The
Particle Identification apparatus is composed of two ring-imaging Cherenkov
detectors that are able to discriminate between different species of charged
hadrons, two calorimeters able to discriminate electrons from hadrons and
muon, and a set of four muon stations downstream the calorimeters.

The size of each event acquired by the LHCb experiment is smaller than for the
other experiments at the LHC [3]. LHCb aims indeed at collecting very large
and pure data samples of $b$ and $c$ hadron decays, to reduce the statistical
uncertainties to a negligible level and to access the rarest decays modes
occurring to less than one hadron in a billion. Both the search for rare
events and the precision measurements on heavy quarks, part of the mainstream
physics programme of the LHCb experiment since its conception, put severe
constraints on the accuracy and precision of the Simulation. For example, the
uncertainties on the efficiencies determined from Simulation should not exceed
the tiny statistical errors obtained with the large samples collected, while
providing a good description of the rejection power of unlikely background
events to effectively design the searches for ultra-rare decays.

As a consequence, very large productions of simulated samples are necessary to
perform several important data analyses. These simulations are performed with
a released and frozen detector configuration and with beam and alignment
conditions tuned to provide the best possible representation of each period of
data taking. It has been estimated that the vast majority of the CPU requests
by the LHCb experiment has been necessary to produce simulated samples, and
that the requests will exceed the pledged resources since 2022 [4].

Important speed-up factors in the simulation can be achieved by reusing the
underlying event associated to the production of a heavy hadron for multiple
decays of the latter [5]. An alternative way forward to reduce the amount of
computing resources needed to produce simulated samples, is to model the
response of each detector element building a sort of cache instead of
simulating the radiation-matter interaction for each particle in each
simulated event.

Machine Learning offers an extremely powerful technique to model the detector
response. In particular, Neural Networks are capable of encoding for each
impinging particle the distribution of possible responses for each detector
element, retaining the capability of modelling the tails in the detector
response giving origin to rare background events, while drastically reducing
the CPU time necessary to simulation.

The rest of this document is organised as it follows. In Section 1, I will
briefly review the concept of Detailed Simulation to establish the baseline on
which injecting the Generative Models for the energy depositions in the
calorimeters and for higher level reconstructed quantities as described in
Sections 2 and 3, respectively. In Sections 4, I will discuss how real data
can be used to train the generative models, obtaining a fast simulation
independent of, and therefore competitive with, the detailed simulation.
Additional details on the training of the Generative Models for detector
simulation are given in Sections 5. Some strategies for the deployment of the
models in the experiment software stack are reviewed in Section 6. Finally, I
report on the comparison between reconstructed data and ultra-fast simulation
in Section 7. Section 8 concludes this paper with a short summary.

##  1 The _Detailed Simulation_ of the LHCb Experiment

Figure 1: Flowchart of the Detailed Simulation processing

Detailed Simulation provides description based on the physical model of the
detector response to the incident particle. As for most HEP experiments, the
first step in the simulation of the LHCb experiment [6] is to simulate the
proton-proton collisions with a Monte Carlo generator such as PYTHIA [7], EPOS
[8], or BcVegPy [9]. The simulation of the heavy hadrons produced in the
collision is obtained running EvtGen [10]. The particles produced in primary
collisions and in subsequent decays are then propagated through the detector
and the energy deposited for radiation-matter interaction is computed with
GEANT4 [11]. The energy deposits (or _hits_) are then digitized [12] and the
hardware trigger is emulated [13, 14] in order to provide to the event
reconstruction pipeline an input as similar as possible to real data. The
reconstruction pipeline is composed of several algorithms dealing with
_tracking_ , associating sets of aligned _hits_ to trajectories of charged
particles (or _track_ s), the _charged particle identification_ , assigning to
each track a likelihood (or other test statistics) for each mass hypothesis,
and _neutral object identification_ mainly intended to identify photons
impinging on the calorimeter. A final step in the reconstruction consists of
algorithms considering sets of tracks to _tag_ some special property of the
event. For example flavour tagging is a technique based on the statistical
study of the particles produced in combination with a $b\bar{b}$ quark pair
developed to discriminate hadrons produced as $B^{0}$ or $\bar{B}^{0}$ and
study flavour oscillations [15, 16, 17, 18].

The reconstructed events are then processed with a centralised selection
framework splitting the events in streams depending on the area of the LHCb
physics programme. The streams are finally processed by the LHCb members
individually to perform statistical analyses on sets of events and draw
conclusions on the underlying physics model.

The procedure to obtain _Detailed Simulation_ samples is depicted as a
flowchart in Figure 1.

##  2 Models for the energy depositions in the calorimeters

The most common application of Machine Learning to Simulation is predicting
the energy deposits of particles in the electromagnetic calorimeter. Indeed,
the detailed simulation of the electromagnetic and hadronic showers in the
calorimeters is computational espensive, accounting for more than half of the
total CPU demand from the LHC experiments [3, 19]. Intense research activity
to speed-up the calorimeter simulation is attested by the many publications on
the subject from the ATLAS [20, 21] and CMS [22, 23, 24] Collaborations, but
also in the context of the effort towards GeantV [25, 26]. A software package
dedicated to the parametric simulation of the calorimeters of ATLAS, CMS and
experiments at future colliders is also being developed under the name of
Delphes [27].

Both traditional parametrization techniques and methods based on Machine
Learning have been explored. When machine learning is employed, it often
relies on Deep Convolutional Generative Adversarial Neural Networks (DCGANs)
[28] to predict the energy deposited in each active volume (voxel) by taking
into account the correlation between adjacent regions. For example, the model
developed for the LHCb experiment solutions based on conditional DCGANs [29]
and a combination of conditional DCGANs and Variational Autoencoders (VAEs)
[30]. In both models the idea is to take as an input the parameters of each
particle (either a charged one or a photon) such as its momentum and slopes,
and translate them to a randomly generated _image_ of the shower in the
electromagnetic calorimeter. The procedure can then be repeated for each track
in the event to reconstruct the overall response of the calorimeter.

These methods provide an output which is intended to be processed with the
whole reconstruction pipeline. The resulting simulated samples contain all the
features exactly as the detailed simulation, including, for example, the
tagging information based on collective features. On the other hand, achieving
good representation of the reconstructed quantities is challenging because the
reconstruction algorithms are not known to the training algorithms, and
therefore the reproduction of reconstructed quantities can be optimised only
indirectly.

##  3 Models for the reconstructed quantities

To achieve a further speed-up of the Simulation improving the quality of the
reconstructed quantities, Musella and Pandolfi proposed to use GANs to
simulate directly the variables obtained from the reconstruction pipeline
[31]. This approach allows for a faster and more accurate simulation with
respect to the simulation of the energy deposits, but the multi-particle
effects on the reconstruction becomes more difficult to take into account
because the reconstruction algorithms are parametrise themselves. As a
consequence, the high-level physics quantities obtained through the _tagging_
part of the reconstruction pipeline represent an additional challenge in this
so called _ultra-fast simulation_.

Nonetheless, the method had success enabling fast simulation of detectors
beyond calorimeters, providing for example the output of the LHCb RICH
detectors [32] and muon system [33]. In both cases, a conditional GAN has been
trained to reproduce the PID-related observables (or, more technically, the
logarithmic likelihood ratio between the muon, kaon, or proton hypothesis and
the pion hypothesis). The reconstruction algorithm used to obtain these
quantities is rather complex, involving a event-wise modelling of the energy
deposits in the calorimeter under various combinations for the mass hypotheses
[34], but the effect of inter-particle correlations on the likelihoods, can be
effectively modelled by including the detector occupancy among the GAN
_conditions_ together with track momentum and pseudorapidity.

##  4 Training on simulated or real data

The most common procedure to train Generative Models to parametrise the
detector response is to generate samples through the Detailed Simulation
pipeline and then train the model to associate generator-level variables to
reconstructed quantities. This simple technique, however, is intrinsically
limited to the quality of the Detailed Simulation or, in other words, a fast
simulation obtained training on detailed simulation will be at most as
accurate as detailed simulation.

As an interesting alternative, the Generative Models can be trained on real
data if following two conditions are satisfied.

  * •

_The acquired dataset is unbiased._ The training data should be selected
without applying selection criteria that, directly or indirectly, modify the
distributions of the reconstructed quantities beyond the effect taken into
account by the _conditions_ passed as inputs to the Generative Model (such as
the kinematics or the detector occupancy). For example, a very common bias is
introduced by trigger algorithms designed to select muons. A training sample
obtained without excluding candidates selected by these lines will result into
Generative Models overestimating the probability of misidentification of
hadrons as muons, because the training sample was enriched of hadrons with a
detector response more similar to that of a muon than the average hadron. In
LHCb this is avoided providing dedicated trigger lines and a dedicated data
processing scheme for the calibration data used to train the Generative Models
for the ultra-fast simulation [35].

  * •

_The generative model supports background subtraction._ Calibration samples
are usually as abundant and pure as possible, but in order to avoid selection
bias, the selections are often loose and some background contribution is
expected. For all the models with a loss function derived more or less
directly from the maximum likelihood principle, it was demonstrated that
weighting the events with properly-computed weights the effect of the
background component on the trained network is null [36, 37]. These
background-subtracting weights, sometimes called ${}_{s}\mathcal{W}$eights in
the context of High Energy Physics, are computed taking into account the
models for the signal and the background components in a discriminant variable
(_e.g._ the invariant mass of the combination of reconstructed particles)
assumed to be independent of all the other (target) variables involved in the
background subtraction [38].

It is worth noticing that an ultra-fast simulation trained on real data is
subject to systematic errors in the description of the detector response which
are completely independent of the detailed simulation. Detailed and ultra-fast
simulation can then cross-check each other, while providing a measure of the
simulation accuracy depending on the specific phase space relevant to each
analysis.

##  5 Automatic hyperparameter optimisation for generative models

The training of GANs is well known within the Machine Learning community to be
a difficult task [39]. Reaching the Nash equilibrium between the Generator and
the Discrminator requires careful tuning of the hyperparameters, such as the
learning rate or the network structure. Unfortunately, the absolute value of
the loss function used to train adversarial models is not meaningful: a
discriminator that finds more difficult to separate real and fake candidates
may mean that the generator is better (and so is the accuracy of the generated
sample) or that the disciminator is worse (and so is the accuracy of the
generated sample).

In order to be able to use Bayesian hyperparameter optimisation frameworks,
such as optuna [40] or scikit-optimize [41], LHCb adopted a _goodness-of-fit_
method based on the ability of a BDT classifier to distinguish training and
generated data entries [42, 43]. The BDT classifier is a third player,
independent of the adversarial training procedure, and is based on a different
classification method with different underlying hypotheses. To provide an
estimate of the uncertainty on quality metric obtained from the BDT, the same
algorithm is trained several times on the same sample with different random
seed. The average and its uncertainty are used as guidance for optimising the
hyperparameters using a distributed version of optuna, named OptunAPI [44].
With this strategy, it is possible to automate the training of GANs on
multiple GPU nodes making the production of a different model per datataking
period accessible.

##  6 Deployment in the software stacks of the experiments

The training and the validation of the generative models relies on Python
frameworks such as _scikit-learn_ [45] for the pre- and post-processing stages
and either _tensorflow_ [46] or _PyTorch_ [47] for the description of the
neural network model.

The pipeline built in Python must then be deployed in the software stack of
the experiments, which is most often written in C++ and released few times a
year. Both $\mathrm{TenorFlow}$ and PyTorch provide C/C++ APIs to deploy
neural network models, however the build system of these large packages often
conflicts with the build system of the experiment’s software. The same applies
for the multi-thread schedulers, designed to meet the high-throughput
requirements of LHC experiments on one side and to speed-up the evaluation of
the networks on the other.

In the context of the CMS and LHCb collaborations, $\mathrm{TenorFlow}$ has
been succesfully integrated within the experiment software111 see CMSSW-DNN
and GaudiTensorFlow on GitLab.  providing evaluations of arbitrarily
complicated networks with short latency. A limitation derived from this
approach is that all the networks within the experiment software must share
the same version of $\mathrm{TenorFlow}$, which makes the maintenance effort
rather significant. An intestersting alternative is to provide an inference
service external to the experiment software, possibly running on the same
machine and interfaced to the data processing software framework via REST APIs
[48]. This approach effectively decouples production code and the internal
architecture of the server which may be based on whatever framework with
whatever thread scheduler and hardware acceleration as long as it is
consistent with the REST interface. Several servers can run concurrently to
rely on multiple frameworks or framework versions. Due to the large latency,
however, this technique is more suitable to process large events with
thousands of tracks than to produce ultra-fast simulation, where the time
spent per event is of the same order as the REST API latency.

To reduce the latency as much as possible, the ATLAS Collaboration has been
developing LWTNN [49], a C++-compiled generic interface for Deep Neural
Networks. LWTNN is composed of an engine compiled together with the experiment
software and released twice a year, and reads configuration files describing
the network architecture and its weights, to reconstruct the model at runtime.
To make the evaluation even faster, LHCb has been developing
RNNGenerator222see RNNGenerator on GitLab to produce C++ code for the Neural
Network and compile it together with the experiment’s software without any
dependence on external frameworks. On the same line, $\mathtt{scikin}$ 333 see
$\mathtt{scikin}$ on GitHub  is a package designed to convert scikit-learn
models in C and distribute them as optimised single-thread shared objects.
More recently, support to simple Keras [50] feed-forward models was added.
Hence, $\mathtt{scikin}$ is much more limited than LWTNN and RNNGenertor in
terms of supported architectures, but eliminates almost any overhead and
enables the distribution of models, including architecture and weights, as
data objects.

The ultra-fast simulation framework of the LHCb Collaboration was originally
designed to rely on GaudiTensorFlow, but an implementation based on
$\mathtt{scikin}$ was found significantly faster. GaudiTensorFlow is now being
phased out and the ultra-fast simulation migrated to $\mathtt{scikin}$.

##  7 Validation

Figure 2:  Comparison of the tracking-related variables for protons obtained
from tagged $\Lambda_{c}^{+}\to pK^{-}\pi^{+}$ decays in calibration data and
through the whole ultra-fast simulation pipeline as implemented in Lamarr.
Figure reproduced from Ref. [51].

Figure 3:  Comparison of the efficiency of two PID requirements with different
thresholds as obtained from tagged $\Lambda_{c}^{+}\to pK^{-}\pi^{+}$ decays
in calibration data from Lamarr. Figure reproduced from Ref. [51].

The fast-simulation validation campaign is performed comparing data from a
calibration sample with the results from Lamarr, the ultra-fast simulation
package under development within LHCb. To validate the generalisation
properties of the models, different decay models are chosen for the training
and the validation of the simulation pipeline. For example, the proton
identification is trained using tagged protons from $\Lambda^{0}\to p\pi^{-}$
decays while the validation relies on tagged protons from $\Lambda_{c}^{+}\to
pK^{-}\pi^{+}$ decays with the charm hadron produced in the semileptonic decay
$\Lambda_{b}^{0}\to\Lambda_{c}^{+}\mu^{-}X$.

Figure 3 reports the comparison between the calibration data and the
simulation for the response of a neural network classifier trained to identify
ghost tracks and for the logarithm of the impact-parameter $\chi^{2}$,
measuring the consistency of the proton track with the primary vertex. The
former is a variable directly obtained as an output of a Generative Model,
while the latter is computed from the simulated properties of the track and
primary vertex, and their (simulated) uncertainties.

Figure 3 shows instead the selection efficiency of a criterion based on the
response of classifier trained to distinguish the particle mass combining the
response of the whole LHCb Particle Identification apparatus: two RICH
detectors, four Muon stations and the electromagnetic calorimeter. In this
case of the simulation, the response of the classifier is obtained directly as
an output of the generative model. The selection efficiency for two different
thresholds are shown as a function of the proton momentum, demonstrating the
capability of the Generative Model to properly model the dependence of the
detector response on the kinematics of the impinging particle. Similar plots
are produced to check good representation of the dependence on the detector
occupancy and other _conditions_. For a deeper discussion of these and
additional comparison plots, see Ref. [51].

##  8 Conclusion

The simulation of the LHCb experiment is an expensive tasks posing severe
pressure on the pledged computing resources. Several solutions to speed-up the
simulation are being explored, among those Machine-Learning Generative Models
are extremely promising. Solutions have been proposed to model both the energy
depositions in the calorimeters and the analysis-level quantities features,
presenting complementary advantages.

The techniques used to train the model on calibration data enables to go
beyond an approximation of the full simulation, obtaining an ultra-fast
simulation potentially competitive with the full simulation and providing
useful cross-checks.

Both the training and the deployment of Generative Models are active research
subjects to keep improving the quality of the simulated features. I discussed
as an example the Bayesian hyperparameter optimisation based on OptunAPI, and
several deployment strategies with different balance between performance and
support to neural network models.

Lamarr, the ultra-fast simulation framework under development in LHCb,
combines several traditional parametrisations and Generative Models to produce
simulated output as similar as possible to reconstructed data. The comparison
of the simulated quantities with calibration samples not used for the training
confirms good generalisation properties.

Upcoming releases of Lamarr will improve the quality of the simulation and
extend the reach of the simulation by adding several reconstructed quantities,
providing a powerful tool within the wide Simulation palette of the LHCb
experiment.

## References

  * [1] LHCb collaboration, _The LHCb Detector at the LHC_ , _JINST_ 3 (2008) S08005.
  * [2] LHCb Collaboration collaboration, _LHCb Detector Performance_ , _Int. J. Mod. Phys. A_ 30 (2015) 1530022 [1412.6352].
  * [3] T. Boccali, _Computing models in high energy physics_ , _Rev. Phys._ 4 (2019) 100034.
  * [4] LHCb collaboration, “LHCb CPU Usage Forecast.” LHCb-FIGURE-2019-18, 2019.
  * [5] D. Müller, M. Clemencic, G. Corti and M. Gersabeck, _ReDecay: A novel approach to speed up the simulation at LHCb_ , _Eur. Phys. J. C_ 78 (2018) 1009 [1810.10362].
  * [6] LHCb collaboration, _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , _J. Phys. Conf. Ser._ 331 (2011) 032047.
  * [7] T. Sjostrand, S. Mrenna and P.Z. Skands, _A Brief Introduction to PYTHIA 8.1_ , _Comput. Phys. Commun._ 178 (2008) 852 [0710.3820].
  * [8] T. Pierog, I. Karpenko, J.M. Katzy, E. Yatsenko and K. Werner, _EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider_ , _Phys. Rev. C_ 92 (2015) 034906 [1306.0121].
  * [9] C.-H. Chang, J.-X. Wang and X.-G. Wu, _BCVEGPY2.0: A Upgrade version of the generator BCVEGPY with an addendum about hadroproduction of the P-wave B(c) states_ , _Comput. Phys. Commun._ 174 (2006) 241 [hep-ph/0504017].
  * [10] D.J. Lange, _The EvtGen particle decay simulation package_ , _Nucl. Instrum. Meth. A_ 462 (2001) 152.
  * [11] GEANT4 collaboration, _GEANT4–a simulation toolkit_ , _Nucl. Instrum. Meth. A_ 506 (2003) 250.
  * [12] G. Corti, M. Cattaneo, P. Charpentier, M. Frank, P. Koppenburg, P. Mato et al., _Software for the LHCb experiment_ , _IEEE Trans. Nucl. Sci._ 53 (2006) 1323.
  * [13] R. Aaij et al., _The LHCb Trigger and its Performance in 2011_ , _JINST_ 8 (2013) P04022 [1211.3055].
  * [14] LHCb collaboration, _Design and performance of the LHCb trigger and full real-time reconstruction in Run 2 of the LHC_ , _JINST_ 14 (2019) P04013 [1812.10790].
  * [15] LHCb collaboration, _New algorithms for identifying the flavour of ${\mathrm{B}}^{0}$ mesons using pions and protons_, _Eur. Phys. J. C_ 77 (2017) 238 [1610.06019].
  * [16] LHCb collaboration, _A new algorithm for identifying the flavour of $B_{s}^{0}$ mesons at LHCb_, _JINST_ 11 (2016) P05010 [1602.07252].
  * [17] LHCb collaboration, _$B$ flavour tagging using charm decays at the LHCb experiment_, _JINST_ 10 (2015) P10005 [1507.07892].
  * [18] LHCb collaboration, _Opposite-side flavour tagging of B mesons at the LHCb experiment_ , _Eur. Phys. J. C_ 72 (2012) 2022 [1202.4979].
  * [19] C. Bozzi, _LHCb Computing Resources: update of 2021 requests, LHCb-PUB-2020-001_, Tech. Rep. (2020).
  * [20] ATLAS collaboration, _AtlFast3: the next generation of fast simulation in ATLAS_ , 2109.02551.
  * [21] C. Krause and D. Shih, _CaloFlow: Fast and Accurate Generation of Calorimeter Showers with Normalizing Flows_ , 2106.05285.
  * [22] M. Paganini, L. de Oliveira and B. Nachman, _CaloGAN : Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks_ , _Phys. Rev. D_ 97 (2018) 014021 [1712.10321].
  * [23] R. Rahmat, R. Kroeger and A. Giammanco, _The fast simulation of the CMS experiment_ , _J. Phys. Conf. Ser._ 396 (2012) 062016.
  * [24] CMS collaboration, _The fast simulation of the CMS detector at LHC_ , _J. Phys. Conf. Ser._ 331 (2011) 032049.
  * [25] G. Amadio et al., _GeantV Alpha Release_ , _J. Phys. Conf. Ser._ 1085 (2018) 032037.
  * [26] G.R. Khattak, S. Vallecorsa, F. Carminati and G.M. Khan, _Fast Simulation of a High Granularity Calorimeter by Generative Adversarial Networks_ , 2109.07388.
  * [27] DELPHES 3 collaboration, _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057 [1307.6346].
  * [28] A. Radford, L. Metz and S. Chintala, _Unsupervised representation learning with deep convolutional generative adversarial networks_ , 1511.06434.
  * [29] V. Chekalina, E. Orlova, F. Ratnikov, D. Ulyanov, A. Ustyuzhanin and E. Zakharov, _Generative Models for Fast Calorimeter Simulation: the LHCb case_ , _EPJ Web Conf._ 214 (2019) 02034 [1812.01319].
  * [30] F. Sergeev, N. Jain, I. Knunyants, G. Kostenkov and E. Trofimova, _Fast simulation of the LHCb electromagnetic calorimeter response using VAEs and GANs_ , _J. Phys. Conf. Ser._ 1740 (2021) 012028.
  * [31] P. Musella and F. Pandolfi, _Fast and Accurate Simulation of Particle Detectors Using Generative Adversarial Networks_ , _Comput. Softw. Big Sci._ 2 (2018) 8 [1805.00850].
  * [32] LHCb collaboration, _Fast Data-Driven Simulation of Cherenkov Detectors Using Generative Adversarial Networks_ , _J. Phys. Conf. Ser._ 1525 (2020) 012097 [1905.11825].
  * [33] G. Sassoli, _Generative Adversarial Networks for the Ultra-Fast Simulation of the Muon Detector of the LHCb experiment at CERN_ , _Thesis discussed in July 2019, University of Florence_ .
  * [34] LHCb RICH Group collaboration, _Performance of the LHCb RICH detector at the LHC_ , _Eur. Phys. J. C_ 73 (2013) 2431 [1211.6759].
  * [35] R. Aaij et al., _Selection and processing of calibration samples to measure the particle identification performance of the LHCb experiment in Run 2_ , _EPJ Tech. Instrum._ 6 (2019) 1 [1803.00824].
  * [36] Y. Xie, _sFit: a method for background subtraction in maximum likelihood fit_ , 0905.0724.
  * [37] M. Borisyak and N. Kazeev, _Machine Learning on data with sPlot background subtraction_ , _JINST_ 14 (2019) P08020 [1905.11719].
  * [38] M. Pivk and F.R. Le Diberder, _SPlot: A Statistical tool to unfold data distributions_ , _Nucl. Instrum. Meth. A_ 555 (2005) 356 [physics/0402083].
  * [39] I.J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair et al., _Generative adversarial nets_ , in _Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2_ , NIPS’14, (Cambridge, MA, USA), p. 2672–2680, MIT Press, 2014.
  * [40] T. Akiba, S. Sano, T. Yanase, T. Ohta and M. Koyama, _Optuna: A next-generation hyperparameter optimization framework_ , 2019.
  * [41] T. Head, M. Kumar, H. Nahrstaedt, G. Louppe and I. Shcherbatyi, _scikit-optimize/scikit-optimize_ , Sept., 2020.  10.5281/zenodo.4014775.
  * [42] M. Williams, _How good are your fits? Unbinned multivariate goodness-of-fit tests in high energy physics_ , _JINST_ 5 (2010) P09004 [1006.3019].
  * [43] C. Weisser and M. Williams, _Machine learning and multivariate goodness of fit_ , 1612.07186.
  * [44] M. Barbetti, _mbarbetti/optunapi_ , _Zenodo_ (2021) 5538989.
  * [45] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel et al., _Scikit-learn: Machine learning in Python_ , _Journal of Machine Learning Research_ 12 (2011) 2825.
  * [46] TensorFlow Developers, _Tensorflow_ , _Zenodo_ (2021) 5189249.
  * [47] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan et al., _Pytorch: An imperative style, high-performance deep learning library_ , _Advances in Neural Information Processing Systems_ (2019) 8024.
  * [48] V. Kuznetsov, L. Giommi and D. Bonacorsi, _MLaaS4HEP: Machine Learning as a Service for HEP_ , _Comput. Softw. Big Sci._ 5 (2021) 17 [2007.14781].
  * [49] D.H. Guest, J.W. Smith, M. Paganini, M. Kagan, M. Lanfermann, A. Krasznahorkay et al., _lwtnn/lwtnn: Version 2.12.1_ , _Zenodo_ (2021) 5082190.
  * [50] F. Chollet et al., “Keras.” https://keras.io, 2015.
  * [51] LHCb collaboration, “Performance of the Lamarr Prototype: the ultra-fast simulation option integrated in the LHCb simulation framework.” LHCb-FIGURE-2019-17, 2019.