Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany<EMAIL_ADDRESS>
# Monolithic MHz-frame rate digital SiPM-IC with sub-100 ps precision and 70
$\mu$m pixel pitch

I. Diehl,Corresponding author  K. Hansen  T. Vanat  G. VignolaAlso at
University of Bonn, Germany  F. Feindt  D. RastorguevAlso at University of
Wuppertal, Germany  and S. Spannagel

###### Abstract

This paper presents the design and characterization of a monolithic integrated
circuit (IC) including digital silicon photomultipliers (dSiPMs) arranged in a
32$~{}\times~{}$32 pixel matrix at 70$~{}\mu$m pitch. The IC provides per-
quadrant time stamping and hit-map readout, and is fabricated in a standard
150-nm CMOS technology. Each dSiPM pixel consists of four single-photon
avalanche diodes (SPADs) sharing a quenching and subsequent processing
circuitry and has a fill factor of 30$~{}\%$. A sub-100 ps precision, 12-bit
time-to-digital converter (TDC) provides timestamps per quadrant with an
acquisition rate of 3 MHz. Together with the hit map, the total sustained data
throughput of the IC amounts to 4 Gbps. Measurements obtained in a dark,
temperature-stable environment as well as by using a pulsed laser environment
show the full dSiPM-IC functionality. The dark-count rate (DCR) as function of
the overvoltage and temperature, the TDC resolution, differential and integral
nonlinearity (DNL/INL) as well as the propagation-delay variations across the
matrix are presented. With aid of additional peripheral test structures, the
main building blocks are characterized and key parameters are presented.

Front-end electronics for detector readout; Particle tracking detectors;
Photon detectors for UV, visible and IR photons (solid-state) (PIN diodes,
APDs, Si-PMTs, G-APDs, CCDs, EBCCDs, EMCCDs, CMOS imagers, etc)

##  1 Introduction

Silicon photomultipliers (SiPMs) are composed of an array of single-photon
avalanche diodes (SPADs), which are photodiodes reverse biased above breakdown
voltage, operating in Geiger-mode. SiPMs working principles and main
advantages like large intrinsic gain (typically 105 to 106), insensitivity
against magnetic fields, and large dynamic range starting with single-photon
response at low bias voltages, are well explained and specified in literature
[1, 2, 3]. Arrays of SPADs are highlighted as ideal candidates when high
sensitivity is required together with high frame rate and precise timing
resolutions [4]. As stated in [5], the capability of photon counting makes
SPADs the detector of choice for applications in which conventional
photodiodes and charge-coupled devices cannot be used, and a large number of
applications for SPAD arrays are listed. Beside the usage as imaging device,
e.g. in [6], also light detection and ranging [7], and direct minimum ionizing
particle detection [8] are mentioned.

In traditional analog SiPMs, each single SPAD is connected to its quenching
resistor forming a so-called microcell. All microcells are connected in
parallel and the common current output signal has to be amplified and
digitized. The number of microcells defines the dynamic range. In digital
SiPMs (dSiPMs) one benefits from the inherent digital behavior of the device,
when sensing directly the output of an individual SPAD by an embedded
quenching and recharging circuitry together with a simple discriminator, e.g.
an inverter. This concept was initially introduced in [9], where SPADs are
integrated with conventional CMOS circuits on the same substrate. This
additional circuitry, in pixel or at the periphery can be used to acquire,
store and transmit data. It negatively impacts the fill factor [3], but offers
the information which SPAD is hit. Individual readout permits the definition
of spatial granularity, e.g. grouping of SPADs in pixel, reduces circuit
complexity and increases the fill factor. Furthermore, it offers temporal
granularity by delivering the timestamp for a section of the pixel array, and
it provides the opportunity to identify and switch off noisy pixels, to adapt
the hold-off time minimizing the afterpulsing probability, and to count hits
within a pixel. The main drawback of implementing SPADs in standard CMOS
technologies compared to a custom technology is the relatively high dark
noise, expressed in terms of dark-count rate (DCR) [10]. To compensate these
DCR and fill-factor issues for the detection of charged particles, one could
implement a dual layer structure of different SPAD arrays for a coincidence
measurement, like in [11]. DCR is one of the most important properties of
SPADs, besides photon-detection probability (PDP), representing the avalanche
probability of the device in response to a photon absorption at a given
wavelength, and afterpulsing, which introduces false events that are
correlated in time to previous detection [3]. Also optical crosstalk, fill
factor, timing resolution and deadtime have to be considered to characterize
SPADs.

The presented dSiPM-IC uses a similar readout concept like in the preceding
project published in [12, 13]. The concept comprises a 32 $\times$ 32 dSiPM-
pixel matrix subdivided into quadrants. Each quadrant consists of a 16
$\times$ 16 dSiPM-pixel array sharing a time-to-digital converter (TDC) and a
validation logic. The 12-bit TDC is designed to provide the timestamp of the
fastest pixel with a time resolution of less than 100 ps. The validation logic
discards undesirable events by setting a threshold. For example, if only one
pixel is hit at the same time, it will be most likely a dark event. Along with
the timing information, also the hit map can be read out continuously in a
frame-based mode at 3 MHz. The previous readout chip was designed in Global
Foundry’s 130-nm CMOS technology comprising a single quadrant. The major
difference between the old and new concept is the interconnection approach of
sensor and readout electronics. In the old concept, the readout IC is flip-
chip connected to the sensor chip utilizing 30 $\mu$m solder spheres at 50
$\mu$m pitch [14]. This hybrid approach enables the direct readout of each
sensor pixel by its corresponding pixel electronics occupying the same pixel
area at about 50 $\mu$m pitch, without affecting the fill factor. Because the
sensitive area of the sensor chip is now vis-á-vis with the readout chip,
direct light illumination of the matrix is not possible. In [13], DCR
measurements were carried out on the first hybrid samples, and the TDC
characterized. In concluding, it was planed to reduce the high DCR by a new
sensor design. A further required step was the redesign of the IC in another
CMOS process. It was decided to follow the current trend and select a
monolithic approach [15]. The chosen process provides on the one hand fully
characterized SPADs in four configurations with suitable properties, and on
the other hand the possibility to redesign the previous developed readout
electronics in a comparable process. The advantage of having sensor and
readout electronics on the same device is the designer’s independency of
sensor provider and interconnection issues. This allowed a fast and cost-
efficient realization of a first proof-of-principle device.

The realized IC is shown and its essential components are described in section
2. Section 3 summarizes the results obtained from measurements and compares
results from test-circuit blocks with results from earlier prototypes [13].
Furthermore, the used Caribou data acquisition (DAQ) system [19] and the test
setups are introduced. The conclusions are given in section 4.

##  2 Chip architecture

###  2.1 Matrix and periphery

The monolithic IC (cf. figure 1a) has been designed in LFoundry’s 150-nm CMOS
technology to take advantage of an available add-on library including single
p+/nwell SPADs with an active area of 20 $\times$ 20 $\mu$m2. This area is
framed by a cathode ring plus an additional pwell ring to minimize the
crosstalk to neighboring cells. The dSiPM-pixel design comprises four SPADs in
parallel sharing the pwell ring and using it partially for the NMOS
transistors of the common pixel electronics to reduce dead area. In this way,
the pixel covers an area of 69.6 $\times$ 76 $\mu$m2 with a fill factor of 30
$\%$. Figure 1 shows the layout of the IC (a) and a photograph of the dSiPM
pixel (b).

[width=0.93]Layout a)b)

Figure 1: a) IC layout (3400 $\times$ 3300 $\mu$m2). b) dSiPM-pixel photograph
(69.6 $\times$ 76 $\mu$m2).

The IC consists of four identical units (quadrants), each one (cf. figure 2)
with a 16 $\times$ 16 dSiPM-pixel matrix, a single 12-bit TDC for event-time
stamping, a validation logic with adjustable settings for discarding
undesirable events, and serializer circuits followed by links (MUX + TX) for
fast (about 1 Gbps) and sustained data readout. In this way, the system
provides the full hit map of a 32 $\times$ 32 dSiPM-pixel matrix together with
four timestamps for each frame.

[width=1]Readout

Figure 2: Quadrant block diagram.

Figure 2 shows a simplified block diagram of one quadrant with all components,
as well as additional blocks used globally for all quadrants. When an
avalanche process is initiated in a SPAD, the pixel draws current from a fast
wired-OR connection, and the earliest pixel in the readout frame triggers the
TDC. Simultaneously, the row-wise wired-ORs are monitored by a validation
logic generating a valid bit for this event. The latched outputs of a
peripheral 40-bit frame counter (FC) complete the quadrant timing information
with the frame number of the event. The TDC and FC are started by the
_Shutter_ signal defining the beginning of a measurement. The 40-bit FC allows
a total recording time of about 100 hours for a measurement. The timing data
are serialized and multiplexed together with the timing data of the other
quadrants. The counted hits in each pixel are serialized row-wise and
multiplexed column-wise to deliver the hit map of the quadrant. In the
periphery, an initialization register provides the control signals for pixel
masking, validation thresholds and TDC settings. Furthermore, clock dividers
and buffers serve to distribute all essential clocks to the components. The
main clocks are provided by the DAQ system. The frame-based readout and
operation mode runs at a targeted frequency of 3 MHz (_FRAME CLK_).
Additionally, a 408-MHz _SYSTEM CLK_ (multiple of the _FRAME CLK_) is required
as reference clock for the TDCs (cf. _TDC Ref Clock_) as well as for the
multiplexers to enable the sustained readout of the entire hit matrix.

###  2.2 Block details and timing

[width=.9]PixelValid a)c)[width=1]TDC b)

Figure 3: Simplified equivalent circuits of: a) the dSiPM pixel, b) the TDC,
and c) the validation logic.

Figure 3a shows the dSiPM-pixel electronics consisting of the four SPADs
connected in parallel, a front-end circuitry using 3.3-V NMOS transistors, and
a readout circuitry operating at core voltage of 1.8 V. The front end allows
an overvoltage (Vov) of maximally 3.6 V, which sets the anode voltage. The
used quenching circuit is similar to the presented one in [16], where the
quenching is performed by a globally biased transistor (cf. _V Quench_) and an
inverter (cf. INV) as comparator for the digital pulse shaping (cf. _Out_). A
clamping transistor (cf. _V Clamp_) limits the inverter input to 1.8 V. The
pixel electronics include the possibility to mask a pixel (cf. _Masking_) via
a single SRAM cell. The hold-off time of the pixel can be adapted by the
global bias voltage _V Quench_. Detected hits can be counted within the
acquisition window by the 2-bit hit counter. If this function is not needed,
the counter can be set as buffer by the control signal _/Set $\\_$2bit_.

Figure 3b shows the TDC, which is divided into a fine and coarse converter.
The fine TDC consists of a delay-locked loop (DLL) with 32 differential delay
elements (DEs). A phase-frequency detector compares the incoming reference
signal (cf. _Rf_) with the delayed one (cf. _Dl_) and steers the following
charge pump like in [17]. The resulting voltage over the loop-filter
capacitance (CLF) together with a calibration circuit controls the delay of
all DEs (cf. _V c_). The DE outputs are latched by the wired-OR trigger signal
until the next frame starts. The latched signals are encoded as thermometer
code and have to be decoded into binary with the aid of a 32-to-5-bit encoder.
After each DLL cycle, a subsequent ripple counter is incremented and serves as
7-bit coarse TDC.

The row-wise wired-OR outputs (cf. _R¡0:15¿_) are monitored by a 4-step
validation logic generating a valid bit for an event, as displayed in figure
3c. In the first validation step, every row output is connected by a
selectable AND/OR gate with its neighboring row output. In each subsequent
step all gates in a column can be programmed as AND or OR (by _Valid
$\\_$cntr¡0:3¿_). In this way, pixel in a cluster of few pixels firing
simultaneously can be identified.

[width=1]Signals

Figure 4: Timing diagram of the frame-based readout.

Figure 4 shows the underlying timing diagram with input clocks, internal
signals and output data for the hit counter set in 1-bit mode. The acquisition
window is defined by the rising edge of _Frame $\\_$rst_ and the falling edge
of _Read_. Afterwards, all collected data are buffered into the serializers
and read out during the following frame. The dynamic range of the TDC with
$2^{7}$/408 MHz = 313.7 ns corresponds to the acquisition window. By using the
hit counter in 2-bit mode, the readout time has to be adapted to enable full
hit-map readout. Accordingly, the acquisition window extends to two frame-
clock cycles.

###  2.3 Test circuits

In the periphery of the matrix, test circuits are integrated for single-block
characterization (cf. figure 1a). These extra blocks comprise a stand-alone
TDC, as well as a chain of receiver (RX) and transmitter (TX) links using the
low-voltage differential signaling (LVDS) standard. The simplified block
schemes are shown in figure 5. The RX consists of a Schmitt-Trigger circuitry
operated at 3.3 V supply voltage and is followed by a voltage-conversation
stage delivering single-ended output. The TX includes an input stage
comprising an edge aligner, voltage converter and buffers to drive the input
of the subsequent current-switching circuitry. To keep the output common mode
stable about 1.2 V, a common-mode feedback (CMFB) steers the output current.
Both circuits are terminated with 100 $\Omega$ on the chip.

[width=1]RXTX

Figure 5: Simplified block diagrams of the RX and TX.

Additionally, sensor-test structures in different configurations are
implemented. For example, a 3 $\times$ 3 test-SPAD array, where the center
SPAD is connected to the same front end as used for the dSiPM pixel, enables
the monitoring of the digital output of a single SPAD.

##  3 Measurements

###  3.1 Experimental setup

For characterization of the IC, the versatile Caribou readout system is used
enabling fast and low-cost implementation of new solid-state detector
prototypes. It offers open access to hardware, firmware [20] and software [21]
speeding up the test setup development. Caribou mainly consists of a system-
on-chip (SoC) evaluation board (Xilinx ZC706) and a control-and-readout (CaR)
interface board. A field-programmable gate array (FPGA) runs custom hardware
blocks for data processing and detector controlling, and an embedded CPU runs
the DAQ and control software. The CaR board provides the physical interface
between the SoC and the detector. It also includes power supplies, analogue
and digital I/Os, a clock generator, analog-to-digital converters, current and
voltage references as well as several connectors. A detector-specific carrier
board comprises the IC, LVDS repeater, and components for filtering,
decoupling and line termination, as well as a trigger logic for an external
pulsed laser source. The carrier board is covered by an aluminum case
protecting the dSiPM-IC from physical damage, external light and acting as a
heat sink. The dark-event based characterization of the dSiPM-pixel matrix was
performed by using a climate chamber and temperature sets between $-$25∘C and
25∘C.

In order to determine the propagation delay along the wired-OR lines, a pulsed
1054-nm laser was operated synchronously to the frame clock. The laser is
movable in all dimensions illuminating sequentially the actual enabled pixel
with a spot diameter of $\sim$0.5 mm.

Another IC test board together with the data-timing generator DTG 5334 (3.35
GHz clock, 0.2 ps step size) and the oscilloscope Teledyne LeCroy SDA-760ZI
were used as second test setup for the characterization of the separate links
and the stand-alone TDC at room temperature.

###  3.2 Results

####  3.2.1 dSiPM-pixel matrix

With the Caribou-test setup, several samples were measured. Initially, the DCR
was monitored as a function of the common SPAD-bias voltage (Vbias) by
accumulating dark hits within 10,000 frames. A hit can only be detected by the
dSiPM pixel (cf. figure 3a) when Vbias exceeds the breakdown voltage (Vbd)
plus the threshold voltage of the inverter. Taking this into account, Vbd
could be estimated by taking the voltage at which hits are detected minus the
threshold voltage of the inverter, which is about 0.6 V. The voltage above Vbd
is defined as the overvoltage (Vov).

[width=.47trim=5 0 30 20,clip]VbdvsTemp a)[width=.47trim=25 5 35
0,clip]DCRvsVov b)[width=.47trim=8 0 35 10,clip]DCRvsTemp c)red black
[width=.47trim=25 0 40 10,clip]DCRdistribution d)red black

Figure 6: a) Vbd versus temperature. Average DCR of all pixels: b) versus Vov
for ten temperature settings, and c) versus temperature at Vov = 1 V and 2 V.
d) Cumulative distribution of all pixel-DCRs in the matrix at three
temperatures and Vov = 2 V.

All results of a representative sample are shown in figure 6. Vbd as a
function of temperature is plotted in figure 6a and the average DCR of all
1024 pixels is plotted in figure 6b. This diagram illustrates the strong
dependency of DCR on temperature and Vov. For example, a decrease of Vov from
2 V to 1 V results in a DCR reduction of about 50 $\%$. Further cooling of the
sample reduces the noise behavior enormously. This is nicely seen in figure 6c
for two different overvoltages. While in 6c the average DCR of all pixels
within the matrix is plotted vs. temperature, figure 6d shows the pixel-DCRs
for $-$25∘C, 0∘C and 25∘C as cumulative distribution (100 $\%$: full matrix of
1024 pixels) for Vov = 2 V. This plot illustrates the wide spread over the
whole matrix, ranging from 600 Hz up to 339.3 kHz for 25∘C and 200 Hz up to
16.2 kHz for $-$25∘C. As a comparison, in [22] a very similar distribution for
an array of 32 $\times$ 32 SPADs at 44.64 $\mu$m pitch using the same sensor
cell and process is shown. The range is here between some ten Hz and $\sim$200
kHz. By using our peripheral 3 $\times$ 3 test-SPAD array we were able to
perform the DCR measurements on single SPADs. The results of three SPADs (550
Hz, 1050 Hz and 15.3 kHz) fit into this distribution plot. If we take the data
from figure 6c (tagged with a ”$\color{red}\downarrow$”) and look for these in
figure 6d, one can say that at least 67 $\%$ of the pixels at $-$25∘C and 73
$\%$ of the pixels at 25∘C have DCR values below the average ones. This
percentage means that only a few very noisy pixels boost the average DCR
value, and this behavior grows with temperature. Disabling 10 $\%$ of the
noisiest pixels leads to an improvement of 40 $\%$ in average DCR.

####  3.2.2 dSiPM-pixel electronics

The previous measurements were taken without using the 2-bit hit-counting
functionality in pixel. This feature allows for the determination of the
deadtime of the pixel (quenching and recharging). Figure 7 shows the deadtime
as function of the global bias voltage (VQuench) of the quenching transistor
(cf. figure 3a), exemplary for a pixel in the center of the laser spot. For
this measurement, two laser pulses were send within a frame starting with
maximal distance (about 450 ns) to each other and minimal VQuench. Below
VQuench = 0.55 V, the electronics counts only one hit, that means the pixel
needs more time for recharging than the acquisition window allows. After first
recognition of the second hit, the deadtime can be defined as distance between
both input pulses. This distance can be decreased with increasing VQuench. The
minimal deadtime is achieved at VQuench = 0.95 V with about 22 ns.

[width=0.8trim=20 0 20 30,clip]Deadtime

Figure 7: Deadtime versus VQuench, with Vov = 2 V at room temperature.

####  3.2.3 Quadrant TDC

[width=0.98]ToA a)b)

Figure 8: a) Histograms of all four quadrant-TDC data of the dark-event
measurements at Vov = 2 V and $-$25∘C, and their Gaussian fit. b) Collected
histograms of all TDC data for $-$25∘ to 25∘C.

In addition to the hit map, the individual quadrant-TDC data are monitored
during the dark-event measurements. In figure 8, time-of-arrival (ToA) is
plotted into histograms measured at Vov = 2 V and different operation
temperatures. At $-$25∘C (figure 8a), the entries for all ToA values are
roughly uniformly distributed over the dynamic range (cf. the Gaussian fit on
the left of figure 8a). With increasing temperature, the entries shift to
lower timestamp values (cf. figure 8b). This is mainly caused by the
temperature-dependent DCR. The probability for more entries at lower
timestamps increases with number of firing pixels within a quadrant, because
all of them share the same TDC, but only the fastest one defines the
timestamp.

[width=1trim=20 30 20 0,clip]fineTDC25C a)b)c)[width=1trim=20 460 20
40,clip]fineTDCm25C d)

Figure 9: a) Histograms, b) DNL, and c) INL of all quadrant fine-TDC data of
the dark-event measurements at Vov = 2 V and 25∘C. d) Histograms at $-$25∘C.

Due to the non-uniform distribution it is not possible to determine the
integral and differential nonlinearity (INL and DNL) of the TDCs over their
total dynamic range with aid of the statistical code density test. Therefore,
figure 9a only shows the histograms of the fine-TDC values. The deviation of
entries for every bin from the mean value of all entries defines the DNL of
each bin. The INL can be calculated as the cumulative sum of DNLs [18]. Figure
9b and c show the DNL and INL for Vov = 2 V and 25∘C, respectively. As visible
in the histogram, the fine TDCs do not tap their full potential of 32 bins (5
bits). That means, the delay of the DEs in the DLL is too large. We believe,
that the circuit is processed in a process corner, where all transistors work
very slowly. Process, voltage and temperature variations can have a big impact
on chip functionality. Corner simulations show a delay range (bin width) of 70
ps until 110 ps. To tune these values two control switches (_TDC
$\\_$cntr¡0:1¿_) are included in the calibration circuit of the TDC (cf.
figure 3b). But for the slow corner, the additional switched current is
insufficient, and only a bin width of about 94 ps is achieved instead of
typical 76.5 ps. Furthermore, the distribution of clocks over the matrix can
provoke run-time discrepancies. This can be the reason for the different
number of fine-TDC bins for the four quadrant TDCs.

For the characterization plots in figure 9, first the DNL standard deviation
(cf. $\sigma$DNL in figure 9b) is determined excluding the last two bins. A
5-$\sigma$ limit on $\sigma$DNL defines the maximal bin number of the fine TDC
or the minimal bin width of the last bin, respectively. Second, the DNL and
INL are determined for the updated fine-TDC bins. This procedure was done for
all dark-event measurements leading to similar results. Here, the effect of
temperature variations is visible in the reached number of bins. The DE delay
decreases with increasing temperature. At $-$25∘C, the dynamic range is
increased by about two bins (cf. figure 9d). This effect makes a stable
temperature environment desirable, and the differences in bin width have to be
considered for the ToA calculations.

[width=1]relFineBinMaps

Figure 10: Offset map (left), and standard deviation of the offset map
(right). Pixel 6/14 (dark blue) was turned off during measurements and its
data are neglected.

As mentioned in section 2, all pixel outputs in a quadrant are connected via
wired-ORs triggering the TDC. These metal interconnections come along with
parasitic elements implicating propagation delays across the matrix. These
delays introduce non-negligible offsets per pixel, which must be considered in
timestamp calculations. Figure 10 shows the offset map (left) and a map with
its standard deviations (right), measured by using the laser source. For this
measurement the pixels are sequentially enabled and illuminated, and the TDC
output per pixel is stored. The offset is calculated by subtracting each pixel
value by the value of the closest pixel to the corresponding TDC (e.g. pixel
1/9 for quadrant 3), which is expected to be the lowest one. In this case, the
offset map illustrates the wired-OR routing to the four TDCs (signal
propagating row-wise from middle to the left and right, and then to the middle
of the edge of each quadrant). The maximum offset in the area of 1.12 $\times$
1.22 mm2 (one quadrant) is (3.45 $\pm$ 0.61) LSB ($\sim$326 ps $\pm$ 86 ps).
The minimum offset is $\leq$ 0 LSB caused by some jitter. The standard
deviation map on the right side of figure 10 highlights the ranges, where the
propagation-delay value falls in a LSB change of the TDCs.

####  3.2.4 Stand-alone TDC

In contrast to the previous subsection, where random dark events were used to
characterize the quadrant fine TDCs, a time-ramp signal was created by driving
the trigger input of the stand-alone TDC with step size of 1 ps over the
entire dynamic range. As reference clock, the 408-MHz system clock was used.
The measured average bin width amounts to $\sim$95 ps. The deviation of each
bin width to the average value defines the DNL, and the deviation to an ideal
line fitted into the step curve defines the INL. Figure 11 exemplarily shows
the measured TDC characteristics of one sample. Displayed are the bin width
and DNL versus code number, and the INL versus time. Additionally, all
measured fine-TDC data are plotted into a histogram, and DNL and INL are
determined with the same method like for the quadrant fine TDCs. The results
for all samples are shown in figure 12. The stand-alone TDCs behave very
similar to the quadrant TDCs. The reached bit resolution is 11.67 bits. The
maximum 12 bit could be reached by decreasing the reference clock to 365 MHz
resulting in an average bin width of 86 ps.

[width=0.98trim=40 40 40 30,clip]TDCout

Figure 11: TDC-bin width, DNL and INL (top to bottom) of one sample at 408 MHz
reference clock.

[width=1trim=20 30 20 30,clip]HistogramsTDC

Figure 12: Histogram, DNL and INL (top to bottom) of four sample’s fine-TDC
data at room temperature.

####  3.2.5 LVDS links

The test RX and TX (cf. figure 5) are connected as input and output for a
buffer with high driving strength. The performance of this RX-TX chain was
measured via eye diagrams. Figure 13 shows the eye diagrams at the typical
data rate of 816 Mbps (top) and at maximum data rate of 1.5 Gbps. At 816 Mbps,
1.3 Gbps and 1.5 Gbps, a bit-error-rate (BER) of ¡10-21, ¡10-15 and ¡10-9 was
achieved, respectively.

[width=0.7]Eyes

Figure 13: Eye diagrams for the RX-TX chain at 816 Mbps (top) and 1.5 Gbps
(bottom).

###  3.3 Summary

The key characteristics are summarized in table 1, and values obtained in our
former design are listed for comparison. In case of the SPAD characteristics,
the DCR per $\mu$m2 could be improved by about one order of magnitude. The
large spread in the noise behavior indicates that a masking of individual SPAD
cells in each pixel could help to improve the DCR without losing all four
SPADs per pixel. A higher fill factor could be achieved by using customized
SPAD cells. The simplification of the pixel electronics shared by four SPADs
led to a reduction of the area by about 80 $\%$. Waiving the in-pixel hit-
counting capability would further reduce the area by 50 $\%$. The TDC
characteristics listed in table 1 are based on the typical system clock of 408
MHz. We obtained very similar characteristics compared to our former design.
However, the TDC power and area requirements increased by about 140 $\%$ and
40 $\%$, respectively. This increase is mainly caused by the higher supply
voltage, different process node and some modifications in the encoder. Taking
the current core area into account, a column-level TDC approach is feasible on
the cost of a 16-times higher power consumption and data throughput. The LVDS
links allow for a 60 $\%$ higher speed level than originally targeted (816
Mbps @ 408 MHz). Therewith, the current design permits a maximum frame rate of
about 4.8 MHz. In contrast to our former design, the area requirement is
increased by 64 $\%$ for TX and RX. The reason is mainly the different process
with other resistor sizes. The conservation of 12-bit resolution entails the
adjustment of the ToA-bin width.

Table 1: Key characteristics of the SPAD, pixel electronics, TDC and LVDS links, determined at room temperature. Parameter | This work | [13]
---|---|---
CMOS node (nm) | 150 | 130
SPAD
Pixel pitch (µm) | 70 | 50
Configuration | 32 $\times$ 32 | 16 $\times$ 16
Fill factor ($\%$) | 30 | 90
mean DCR (Hz/($\mu$m2) @ Vov = 1 V | 8.7 | 80
Pixel electronics
Area ($\mu$m $\times~{}\mu$m) | 70 $\times$ 5 + 3 $\times$ 17111cf. figure 1b | 40 $\times$ 45
Power ($\mu$W) | 10 | 25
TDC
Resolution (bit) | 11.67 | 12
Precision (ps) | 95.8 $\pm~{}$13.65 | 77.19 $\pm$ 7.53
max. DNL (LSB) |  $-$0.74 / 0.35 |  $-$0.46 / 0.64
max. INL (LSB) |  $-$1.43 / 1.39 |  $-$1.33 / 0.93
Power (mW) | 11 | 4.6
Area ($\mu$m $\times~{}\mu$m) | 78 $\times$ 157 | 55 $\times$ 160
Links (RX-TX chain)
| @BER=10-15 | 1.3 | 1.2
max. data throughput (Gbps) | @BER=10${}^{-9~{}}$ | 1.5 | 1.6
Power (mW) | RX / TX | 3 / 37 | 8 / 48
| RX | 39 $\times$ 38 | 30 $\times$ 30
Area ($\mu$m $\times$ $\mu$m) | TX | 73 $\times$ 94 | 60 $\times$ 70

##  4 Conclusions

We presented a monolithic 32 $\times$ 32 dSiPM-pixel matrix IC with 70 $\mu$m
pitch and 30 $\%$ fill factor designed and fabricated in LF’s 150-nm CMOS
technology. It enables full hit-map readout and provides sub-100 ps time
stamping for each quadrant. Main focus has been taken on dark-event
measurements in a climate chamber to show the dependency of DCR on temperature
and overvoltage. Also the event-related timestamps for different temperatures
have been analyzed. With aid of a pulsed laser, propagation delays across the
matrix have been determined. Furthermore, the TDC resolution, DNL and INL as
well as data transmission speed limits have been identified and compared with
a previous prototype designed for a hybrid concept.

The authors would like to thank A. Venzmer, E. Wüstenhagen and D. Gorski for
test-board and mechanical case design, test setup, as well as chip assembly.
We are grateful to C. Reckleben and S. Lachnit for fruitful discussions and
manuscript reading.

## References

  * [1] A.N. Otte, B. Dolgoshein, J. Hose, S. Klemin, E. Lorenz, G. Lutz et al., _Prospects of using silicon photomultipliers for the astroparticle physics experiments EUSO and MAGIC_ , https://doi.org/10.1109/TNS.2006.870575 _IEEE Transactions on Nuclear Science_ 53 (2006) 636–640.
  * [2] A.N. Otte, _The silicon photomultiplier - A new device for high energy physics, astroparticle physics, industrial and medical applications_ , in proceedings of the _SNIC Symposium_ , Stanford, CA, U.S.A., 3–6 April 2006.
  * [3] C. Bruschini, H. Homulle, I.M. Antolovic, S. Burri and E. Charbon, _Single-photon SPAD imagers in biophotonics: Review and outlook_ , https://doi.org/10.48550/arXiv.1903.07351 _arXiv:1903.07351_ [physics.ins-det] (2019).
  * [4] F. Guerrieri, S. Bellisai, A. Tosi, G. Padovini, F. Zappa and S. Tisa, _SPAD arrays for parallel photon counting and timing_ , in proceedings of the https://doi.org/10.1109/PHOTONICS.2010.5698906 _2010 23rd Annual Meeting of the IEEE Photonics Society_ , Denver, CO, U.S.A., 7-11 Nov. 2010, 355–356.
  * [5] M.-L. Wu, E. Ripiccini, E. Kizilkan, F. Gramuglia, P. Keshavarzian, C.A. Fenoglio et al., _Radiation hardness study of single-photon avalanche diode for space and high energy physics applications_ , https://doi.org/10.3390/s22082919 _Sensors_ 22 (8) (2022).
  * [6] K. Morimoto, A. Ardelean, M.-L. Wu, A.C. Ulku, I.M. Antolovic, C. Bruschini et al., _A megapixel time-gated SPAD image sensor for 2D and 3D imaging applications_ , https://doi.org/10.48550/arXiv.1912.12910 _arXiv:1912.12910_ [eess.IV] (2019).
  * [7] K. Yoshioka, H. Kubota, T. Fukushima, S. Kondo, T.T. Ta, H. Okuni et al., _A 20-ch TDC/ADC hybrid architecture LiDAR SoC for 240 $\times$ 96 pixel 200-m range imaging with smart accumulation technique and residue quantizing SAR ADC_, https://doi.org/10.1109/JSSC.2018.2868315 _IEEE Journal of Solid-State Circuits_ 53 (2018) 3026–3038.
  * [8] F. Gramuglia, E. Ripiccini, C.A. Fenoglio, M.-L. Wu, L. Paolozzi, C. Bruschini et al., _Sub-10 ps minimum ionizing particle detection with geiger-mode APDs_ , https://doi.org/10.48550/arXiv.2111.09998 _arXiv:2111.09998_ [physics.ins-det] (2021).
  * [9] T. Frach, G. Prescher, C. Degenhardt, R. de Gruyter, A. Schmitz and R. Ballizany, _The digital silicon photomultiplier – principle of operation and intrinsic detector performance_ , in proceedings of the https://doi.org/10.1109/NSSMIC.2009.5402143 _2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC)_ , Orlando, FL, U.S.A., 2009, 1959–1965.
  * [10] G. Torilla, J. Minga, P. Brogi, G. Collazuol, G.-F. Dalla Betta, P. Marrocchesi et al., _DCR and crosstalk characterization of a bi-layered 24 $\times$ 72 CMOS SPAD array for charged particle detection_, https://doi.org/10.1016/j.nima.2022.167693 _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 1046 (2023) 167693.
  * [11] L. Ratti, P. Brogi, G. Collazuol, G.-F. Dalla Betta, P.S. Marrocchesi et al., _Layered CMOS SPADs for Low Noise Detection of Charged Particles_ , https://doi.org/10.3389/fphy.2020.607319 _Front. Phys._ 8, 607319 (2020).
  * [12] I. Diehl, K. Hansen, K. Krüger, C. Reckleben, F. Sefkow, L. Andricek et al., _Readout ASIC for fast digital imaging using SiPM sensors: Concept study_ , in proceedings of the https://doi.org/10.1109/NSSMIC.2015.7581816 _2015 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_ , San Diego, CA, U.S.A., 2015, 1–3.
  * [13] I. Diehl, K. Hansen, K. Krüger, C. Reckleben, F. Sefkow, Y. Yao et al., _Readout of digital SiPMs_ , in proceedings of the https://doi.org/10.1109/NSSMIC.2018.8824395 _2018 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_ , Sydney, NSW, Australia, 2018, 1–3.
  * [14] S. Kousar, K. Hansen, and T. F. Keller, _Laser-Assisted Micro-Solder Bumping for Copper and Nickel-Gold Pad Finish_ , https://doi.org/10.3390/ma15207349 _Materials_ 15 (20) (2022) 7349.
  * [15] J. Sonnefeld, _Design and Performance of HV CMOS Sensors for Future Colliders by the RD50 Collaboration_ , https://doi.org/10.48550/arXiv.2307.08600 _arXiv:2307.08600_ [physics.ins-det] (2023).
  * [16] H. Xu, L. H. C. Braga, D. Stoppa and L. Pancheri, _Characterization of Single-Photon Avalanche Diode arrays in 150nm CMOS technology_ , in proceedings of the https://doi.org/10.1109/AISEM.2015.7066818 _2015 XVIII AISEM Annual Conference_ , Trento, Italy, 2015, 1–4.
  * [17] M. Estebsari,M. Gholami, M. J. Ghahramanpour, _A wide range delay locked loop for low power and low jitter applications_ , https://doi.org/10.1002/cta.2401 _International Journal of Circuit Theory and Applications_ 46 (3) (2018) 401–414.
  * [18] F. Villa, B. Markovic, S. Bellisai, D. Bronzi, A. Tosi, F. Zappa et al., _SPAD Smart Pixel for Time-of-Flight and Time-Correlated Single-Photon Counting Measurements_ , https://doi.org/10.1109/JPHOT.2012.2198459 _IEEE Photonics Journal_ 4 (3) (2012) 795–804.
  * [19] T. Vanat, _Caribou – A versatile data acquisition system_ , in proceedings of https://doi.org/10.22323/1.370.0100 _PoS TWEPP2019_ , Santiago de Compostela, GAL, Spain, 2019, 1–5.
  * [20] _The Caribou DAQ system_ , https://gitlab.cern.ch/Caribouhttps://gitlab.cern.ch/Caribou.
  * [21] _Peary Caribou DAQ framework_ , https://gitlab.cern.ch/Caribou/pearyhttps://gitlab.cern.ch/Caribou/peary.
  * [22] M. Zarghami, L. Gasparini, L. Parmesan, M. Moreno-Garcia,A. Stefanov, B. Bessire et al., _A 32 $\times$ 32-Pixel CMOS Imager for Quantum Optics With Per-SPAD TDC, 19.48$\%$ Fill-Factor in a 44.64-$\mu$m Pitch Reaching 1-MHz Observation Rate_, https://doi.org/10.1109/JSSC.2020.3005756 _IEEE Journal of Solid-State Circuits_ 55 (10) (2020) 2819–2830.


# Measuring $hhWW$ Coupling at Lepton Colliders

Qing-Hong Cao<EMAIL_ADDRESS>School of Physics, Peking University,
Beijing 100871, China  Center for High Energy Physics, Peking University,
Beijing 100871, China  Kun Cheng<EMAIL_ADDRESS>School of Physics,
Peking University, Beijing 100871, China  Yandong Liu<EMAIL_ADDRESS>Key
Laboratory of Beam Technology of Ministry of Education, College of Nuclear
Science and Technology, Beijing Normal University, Beijing 100875, China
Beijing Radiation Center, Beijing 100875, China  Xiao-Rui Wang
<EMAIL_ADDRESS>School of Physics, Peking University, Beijing 100871,
China

###### Abstract

The quartic Higgs-gauge-boson coupling $g_{hhWW}$ is sensitive to the
electroweak symmetry breaking mechanism, however, it is challenging to be
measured at the large hadron collider. We show that the coupling can be well
probed at future lepton colliders through di-Higgs boson production via the
vector boson fusion channel. We perform a detailed simulation of
$\ell^{+}\ell^{-}\to\nu\bar{\nu}hh\to 4b+\not{E}_{T}$ with parton showering
effects at $e^{+}e^{-}$, $\mu^{+}\mu^{-}$ and $e^{-}\mu^{+}$ colliders. We
find that the regions of $g_{hhWW}/g_{hhWW}^{\rm SM}<0.86$ and
$g_{hhWW}/g_{hhWW}^{\rm SM}>1.32$ can be discovered at the $5\sigma$
confidence level at the 10 TeV $\mu^{+}\mu^{-}$collider with an integrated
luminosity of 3 ab-1.

##  I Introduction

The discovery of the Higgs boson Aad  _et al._ (2012); Chatrchyan  _et al._
(2012) has significant implications for the study of its characteristics,
including its mass, spin, and interactions with other particles in the
Standard Model (SM). Yet, our understanding of the electroweak symmetry
breaking (EWSB) mechanisms and the true character of the Higgs boson remains
shrouded in mystery. The couplings of the Higgs boson to $W^{\pm}$ and $Z$
bosons, $g_{hVV}$ and $g_{hhVV}$, pivotal for insights into EWSB, are
intrinsically linked to the interactions of the gauge interactions and pattern
of EWSB. The SM values of the couplings, $g^{\text{SM}}_{hVV}$ and
$g^{\text{SM}}_{hhVV}$, are uniquely fixed by the weak charge of the scalar
doublet. The interrelationship of these couplings manifests the internal
consistency of the SM and ensures perturbative unitarity of gauge boson
scatterings.

These couplings, however, can be altered in new physics (NP) models. To
describe the NP effects in these couplings, we define $\kappa_{hhVV}$ and
$\kappa_{hVV}$ as:

$\displaystyle\kappa_{hhVV}\equiv\frac{g_{hhVV}}{g^{\text{SM}}_{hhVV}},\quad\kappa_{hVV}\equiv\frac{g_{hVV}}{g^{\text{SM}}_{hVV}}.$
(1)

Any deviation of $\kappa_{hhVV}$ and $\kappa_{hVV}$ from unity implies the
existence of NP. For instance, a singlet scalar extension model predicts that
both the $g_{hhVV}$ and $g_{hVV}$ couplings will decrease compared to the SM
due to mixings, leading to $\kappa_{hhVV}$ and $\kappa_{hVV}$ that are less
than one. Doublet scalar extension models might decrease $g_{hVV}$, but leave
$g_{hhVV}$ unchanged due to EWSB constraints. Consequently, in these models,
$\kappa_{hVV}$ is less than one, but $\kappa_{hhVV}$ retains the SM value.
Furthermore, a larger $\kappa_{hVV/hhVV}$ is predicted in the scalar triplet
or larger representations, owing to increasing Clebsch-Gordan (CG)
coefficients. On the other hand, both the couplings $g_{hVV}$ and $g_{hhVV}$
are reduced in the Composite Higgs Model (CHM), in which the coupling ratio
$\kappa_{hVV}=\sqrt{1-\xi}$ and $\kappa_{hhVV}=1-2\xi$, where
$\xi=v^{2}/f^{2}$, $v$ is the electroweak symmetry breaking scale and $f$
represents the compositeness scale Agashe  _et al._ (2005).

Measuring $\kappa_{hhVV}$ and $\kappa_{hVV}$ and investigating their relation
would shed lights on the mechanism of EWSB. The ratio $\kappa_{hVV}$ can be
precisely pinned down by the vector boson fusion (VBF) type single Higgs
production at the LHC Sirunyan  _et al._ (2019); Sharma and Shivaji (2022) as
$\kappa_{hVV}$ can be extracted directly from
$\kappa_{hVV}=\sqrt{\sigma_{h}/\sigma_{h}^{\rm SM}}$, where $\sigma_{h}$ is
the cross section of VBF-type single Higgs production and $\sigma_{h}^{\rm
SM}$ is the SM prediction. However, it is difficult to measure the ratio
$\kappa_{hhVV}$ at the LHC, because $\kappa_{hhVV}$ can only be measured from
VBF type di-Higgs production processes while the gluon fusion processes is the
dominant di-Higgs production processes at the LHC de Florian and Mazzitelli
(2013). Presently, the LHC provides a loose constraint on $\kappa_{hhVV}$ with
$0<\kappa_{hhVV}<2.1$ Aad  _et al._ (2020, 2023).

In our study, we investigate the VBF-type di-Higgs production at future high-
energy lepton colliders Roloff  _et al._ (2020); Han  _et al._ (2021a), and we
focus on the measurement of Higgs-gauge-boson couplings. The rest of this
paper is organized as follows. The cross section of different di-Higgs
production is described in Sec. II. The selection and reconstruction of events
is presented in Sec. III, and the result is shown in Sec. IV. Finally we
conclude in Sec. V.

##  II Di-Higgs production cross section at future lepton colliders

Figure 1: The cross sections of di-Higgs boson production at lepton colliders
as a function of colliding energies ($\sqrt{s}$), where $W$BF and $Z$BF stands
for the VBF process mediated by the $W$ and $Z$ boson, respectively.

At high-energy lepton colliders, di-Higgs boson production can be instigated
through several processes, including the VBF process, the $Z$-boson associated
production ($Zhh$) process, and the loop-induced $hh$ pair production process.
Figure 1 plots the cross section for each process as a function of the
collision energy $\sqrt{s}$. The contribution of loop-induced process is very
small due to the loop suppression, and the $s$-channel production cross
section of $Zhh$ decreases with collision energy when $\sqrt{s}$ exceeds 1 TeV
due to $s$-channel suppression. Contrarily, the VBF production cross section
increases with $\log\sqrt{s}/m_{V}$, as the the collinear splitting mechanism
of the electroweak gauge bosons becomes the dominant phenomena of the initial
state radiation at high energy Han  _et al._ (2021b). Besides, the cross
section of the charged current process (mediated by $W^{\pm}$) is
approximately an order of magnitude larger than that of the neutral current
process (mediated by $Z$). Therefore, we mainly focus on the $W^{\pm}$ fusion
process in the subsequent analysis.

Figure 2: Feynman diagrams of the di-Higgs boson production via VBF processes
at lepton colliders. Figure 3: Contours of $\kappa_{hhWW}$ in the plane of
$\sigma_{h}/\sigma_{h}^{\rm SM}$ and $\sigma_{hh}/\sigma_{hh}^{\rm SM}$, where
$\sigma_{h}$ and $\sigma_{hh}$ presents the cross section of the VBF type
single-Higgs boson production and the di-Higgs boson production at the 10 TeV
muon collider, respectively.

The VBF type di-Higgs production involves both the couplings $g_{hWW}$ and
$g_{hhWW}$; see Fig. 2. Therefore, measuring $g_{hhWW}$ via di-Higgs
production also relies on the value of $g_{hWW}$. Fortunately, the ratio
$\kappa_{hWW}$ can be precisely measured from single Higgs production via VBF
process and $Vh$ associated process. Combining both the single and di-Higgs
production, one can determine the $\kappa_{hhWW}$. In Fig. 3, we plot the
contours of $\kappa_{hhWW}$ in the plane of $\sigma_{h}/\sigma_{h}^{\rm SM}$
and $\sigma_{hh}/\sigma_{hh}^{\rm SM}$. Even though a degeneracy of
$\kappa_{hhWW}$ appears for given total cross sections $\sigma_{h}$ and
$\sigma_{hh}$, it can be resolved by analyzing the $m_{hh}$ distribution. In
the work we performed simulations at three kinds of lepton colliders

  * •

$e^{+}e^{-}$ collider with collision energy of $\sqrt{s}=1$ TeV;

  * •

$e\mu$ collider with collision energies of $\sqrt{s}=1$ TeV and $\sqrt{s}=2$
TeV Lu  _et al._ (2021);

  * •

$\mu^{+}\mu^{-}$collider with collision energy of $\sqrt{s}=$10 TeV Accettura
_et al._ (2023).

Assuming the SM value of trilinear Higgs boson self-coupling, the coupling
dependence of the cross section can be parametrized as

$\displaystyle\sigma^{hh}=$
$\displaystyle\kappa_{hhWW}^{2}\sigma_{a}+\kappa_{hWW}^{4}\sigma_{b}+\kappa_{hWW}^{2}\sigma_{c}+\kappa_{hhWW}\kappa_{hWW}^{2}\sigma_{ab}+\kappa_{hWW}^{3}\sigma_{bc}+\kappa_{hhWW}\kappa_{hWW}\sigma_{ac},$
(2)

where $\sigma_{a}$, $\sigma_{b}$ and $\sigma_{c}$ parametrize the contribution
from the Feynman diagram (a), (b) and (c) in Fig. 2, respectively;
$\sigma_{ab}$, $\sigma_{ac}$ and $\sigma_{bc}$ denote the interference between
each Feynman diagram. For the process of
$l^{-}l^{+}\to\nu_{l}\bar{\nu}_{l}hh$, the contribution from each Feynman
diagram and their interference effects are shown in Fig. 4, and we also show
the values with the benchmark collision energies in Table 1.

Figure 4: Cross section of the $W$-boson fusion process $\ell^{+}\ell^{-}\to\bar{\nu}_{\ell}\nu_{\ell}hh$ from the contribution of each diagram and their interferences.  Table 1: Parametrization of parton level cross sections of $l^{-}l^{+}\to\nu_{l}\bar{\nu}_{l}hh$ at lepton colliders. cross section [fb] | $\sigma_{a}$ | $\sigma_{b}$ | $\sigma_{c}$ | $\sigma_{ab}$ | $\sigma_{bc}$ | $\sigma_{ac}$
---|---|---|---|---|---|---
$1\operatorname{TeV}$ $e^{+}e^{-}$ collider | $0.54$ | $1.41$ | $0.09$ | $-1.69$ | $-0.62$ | $0.35$
$1\operatorname{TeV}$ $e\mu$ collider | $0.54$ | $1.41$ | $0.09$ | $-1.69$ | $-0.62$ | $0.35$
$2\operatorname{TeV}$ $e\mu$ collider | $4.80$ | $9.02$ | $0.34$ | $-12.72$ | $-2.53$ | $1.53$
$10\operatorname{TeV}$ $\mu^{+}\mu^{-}$ collider | $177.2$ | $206.1$ | $0.88$ | $-378.10$ | $-12.11$ | $9.08$

When the ratio $\kappa_{hhWW}$ has a deviation from its SM value with
$\delta\kappa_{hhWW}$, namely $\kappa_{hhWW}=1+\delta\kappa_{hhWW}$, the
leading contribution of the deviation arises from diagram (a) and its
interference between other diagrams; see Eq. (2). It is observed in Fig. 4
that $\sigma_{a}$, $\sigma_{b}$ and $\sigma_{ab}$ increase with collision
energy and dominate at high energy. Therefore, the cross section becomes more
sensitive to $\kappa_{hhWW}$ with a larger collision energy. The underlying
physics of this explosion of $\sigma_{a}$, $\sigma_{b}$ and $\sigma_{ab}$ with
respect to the energy is the violation of perturbative unitarity; see Appendix
A for details. Therefore, the high energy lepton colliders are ideal to test
NP models that change the Higgs-gauge-boson coupling $g_{hhWW}$.

##  III Collider Simulation

We consider three lepton colliders, the $e^{+}e^{-}$ collider and
$\mu^{+}\mu^{-}$collider with symmetric incident beam energies and the $e\mu$
collider with asymmetric incident beams, and we explore their potential of
measuring the $\kappa_{hhWW}$ coupling through the VBF-type di-Higgs
production. Owing to the clean environment of the lepton colliders, one can
reconstruct the two Higgs bosons in the $b\bar{b}$ decay mode, the mode with
the largest branching ratio. This leads to a distinct signature comprising
four $b$-jets and missing energy $\not\mathrel{E}_{T}$ from two invisible
neutrinos, in all the three lepton colliders.

The primary background processes for the $4b+$$\not\mathrel{E}_{T}$ signal are
$\ell_{i}^{-}\ell_{j}^{+}\to\nu_{i}\bar{\nu}_{j}b\bar{b}b\bar{b}$ from on-
shell $hZ/ZZ$ production, and $\ell^{\pm}\gamma\to
b\bar{b}jj\nu\ell(\bar{\nu}\ell)$ from on-shell $hW^{\pm}/ZW^{\pm}$ production
Roloff  _et al._ (2020). In addition, the $\ell^{+}\ell^{-}\to
b\bar{b}b\bar{b}$ process, arising from on-shell $Zh$ or $ZZ$ production, may
also contribute to the background at the $e^{+}e^{-}$ collider and
$\mu^{+}\mu^{-}$collider, where $\not\mathrel{E}_{T}$ is mimicked by
undetected particles in the showering. Nevertheless, the reducible $4b$
backgrounds can be effectively suppressed by employing a relatively large
$\not\mathrel{E}_{T}$ cut or a recoil mass cut, as discussed below.

###  III.1 Symmetric $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$ Colliders

First, we focus on the symmetric colliders, namely the $e^{+}e^{-}$ and
$\mu^{+}\mu^{-}$ colliders. The $e^{+}e^{-}$ collider is proposed as a
precision machine that benefits from the clear collision environment ILC
(2013); Asner  _et al._ (2013); Dong  _et al._ (2018); Aryshev  _et al._
(2022). The $\mu^{+}\mu^{-}$collider, as proposed in Budker (1969);
Parkhomchuk and Skrinsky (1983); Neuffer (1983), offers not only a clean
environment but also a high collision energy Palmer  _et al._ (1996);
Stratakis  _et al._ (2022); Abbott  _et al._ (2022); Aimè  _et al._ (2022);
Han  _et al._ (2022a, b). While technical difficulties Delahaye  _et al._
(2019); Long  _et al._ (2021) arise from the short lifetime, they can be
compensated by the advantage of increasing instantaneous luminosity with beam
energy Long  _et al._ (2021). We use $e^{+}e^{-}$ collider with $\sqrt{s}=1$
TeV and $\mu^{+}\mu^{-}$collider with $\sqrt{s}=10$ TeV as examples in our
analysis.

In the following simulation, we generate parton-level events using MadEvent
Alwall  _et al._ (2014), then pass them to Pythia8 Sjostrand  _et al._ (2008)
for showering and hadronization. The showered events are clustered using
Fastjet Cacciari  _et al._ (2012) using anti-$k_{T}$ algorithm with clustering
parameter $\Delta R=0.5$, and we subsequently simulate collider effects with
Delphes de Favereau  _et al._ (2014). We adopt the tracking and energy
resolution modeling used in Delphes as specified in Dong  _et al._ (2018). We
select a benchmark $b$-jet (mis-)tagging efficiencies from Ref. Dong  _et al._
(2018) with the higher identification rate,

$\displaystyle p_{b\to b}=0.9,~{}~{}~{}p_{c\to b}=0.3,~{}~{}~{}p_{j\to
b}=0.05,$ (3)

where $j$ represents gluon or light quarks.

For the pre-selection, we require both the signal and backgrounds to satisfy
the following cuts:

$\displaystyle
n^{\ell}(p_{T}>10\operatorname{GeV})=0,~{}\not\mathrel{E}_{T}>10\operatorname{GeV},$
$\displaystyle
p_{T}^{jet}>15~{}\operatorname{GeV},~{}-4.0<\eta^{jet}<4.0,~{}\Delta
R^{mn}>0.5,$ (4)

where $\ell$ denotes an electron or muon, $m(n)$ represents a reconstructed
lepton or jet, and $\Delta R$ is the rapidity and azimuthal distance defined
as $\Delta R^{mn}=\sqrt{\Delta\eta_{mn}^{2}+\Delta\phi_{mn}^{2}}$. To
eliminate backgrounds mediated by neutral current interactions, we require a
$\not\mathrel{E}_{T}$ signal and no isolated leptons with transverse momenta
larger than $10\operatorname{GeV}$. To effectively reconstruct the two Higgs
bosons, we require at least four $b$-jets and divide the four leading ones
into two pairs by minimizing
$\chi^{2}_{ijkl}=(m_{b_{i}b_{j}}-m_{h})^{2}+(m_{b_{k}b_{l}}-m_{h})^{2}$. The
two pairs of $b$-jets are treated as two Higgs boson candidates. Then we
impose the requirement that the invariant mass of the two $b$-jets in each
pair must satisfy

$100\operatorname{GeV}<m_{bb}<150\operatorname{GeV},$ (5)

which enhances signal-to-backgrounds ratio dramatically; see Tables 2 and 3.

To further suppress electroweak resonance backgrounds, we require the
likelihood of the four $b$-jets originating from the Higgs boson pair to be
larger than that from a $ZZ$ pair or $Zh$ pair. We achieve this by employing
the function $\chi^{2}(m_{1},m_{2})$, defined as

$\displaystyle\chi^{2}(m_{1},m_{2})\equiv\min_{{i,j,k,l}}\left[(m_{b_{i}b_{j}}-m_{1})^{2}+(m_{b_{k}b_{l}}-m_{2})^{2}\right],$
(6)

which represents the minimal mass difference between two $b$-jet pairs and
$m_{1(2)}$. We then require the four leading $b$-jets in each event to satisfy

$\displaystyle{\rm{HHCUT}:}~{}~{}~{}\chi^{2}(m_{h},m_{h})<\chi^{2}(m_{Z},m_{Z}),~{}~{}~{}\chi^{2}(m_{h},m_{h})<\chi^{2}(m_{h},m_{Z}).$
(7)

The selection criteria employed in this study effectively reduces the
$ZZ$+$\not\mathrel{E}_{T}$ and $Zh$+$\not\mathrel{E}_{T}$ backgrounds and
preserves the signal events, as illustrated in the fourth columns of Tables 2
and 3.

At last, the process $e^{+}e^{-}\to b\bar{b}b\bar{b}$ stemming from $Zh/ZZ$
resonant production often presents a non-negligible background. For example,
at the $1\operatorname{TeV}$ $e^{+}e^{-}$ collider, the cross section of these
two non-intrinsic backgrounds are still larger than the signal cross section
after imposing all the previous cuts. Although it is feasible to reject these
backgrounds by increasing the $\not\mathrel{E}_{T}$ cut in Eq. (4), this
approach also reduces signal events considerably. A more effective strategy is
to impose a cut on the recoil mass, defined as

$M_{\rm recoil}\equiv\sqrt{(p_{1}+p_{2}-p_{h_{1}}-p_{h_{2}})^{2}},$ (8)

where $p_{1}$ and $p_{2}$ represent the momenta of initial lepton, and
$p_{h_{1}}$ and $p_{h_{2}}$ denote the momenta of reconstructed Higgs boson.
We require a hard cut on $M_{\rm recoil}$

$\displaystyle M_{\rm recoil}\geq 200\operatorname{GeV},$ (9)

as the two neutrinos in the signal tend to fly back-to-back and lead to a
large recoil mass. We find that, after the $M_{\rm recoil}$ cut, the non-
intrinsic backgrounds are negligible, and the signal events almost remain the
same; see Table 2. With an integrated luminosity of $3\operatorname{ab}^{-1}$,
the signal strength of the total cross section of di-Higgs production,
$\mu_{hh}=\sigma_{hh}/\sigma_{hh}^{\rm SM}$, can be constrained to
$0.35<\mu_{hh}<1.80$ at the $1\operatorname{TeV}$ $e^{+}e^{-}$ collider and
$0.9<\mu_{hh}<1.1$ at the $10\operatorname{TeV}$ $\mu^{+}\mu^{-}$collider, at
95% confidence level.

Table 2: Cut flow for the SM $\nu\bar{\nu}hh$ di-Higgs signal and backgrounds
cross sections at $1\operatorname{TeV}$ electron-positron collider (left) and
10 TeV muon collider (right).

$\sigma$ [ab] | pre-cuts |  $m_{bb}$ cut | HHCUT |  $M_{\rm recoil}$ cut
---|---|---|---|---
Sig. | 15.9 | 9.7 | 8.3 | 5.7
$\nu\nu hZ$ | 54.4 | 11.1 | 6.3 | 5.5
$\nu\nu ZZ$ | 73.7 | 9.1 | 3.1 | 2.9
$\nu W^{\pm}h$ | 45.2 | 3.9 | 2.3 | 2.3
$\nu W^{\pm}Z$ | 47.0 | 2.5 | 1.1 | 1.1
Bkg. | - | - | - | 11.8

$\sigma$ [ab] | pre-cuts |  $m_{bb}$ cut | HHCUT |  $M_{\rm recoil}$ cut
---|---|---|---|---
Sig. | 484.4 | 261.5 | 226.5 | 226.3
$\nu\nu hZ$ | 1163.1 | 168.6 | 97.8 | 97.6
$\nu\nu ZZ$ | 1557.9 | 121.2 | 36.7 | 36.6
$\nu W^{\pm}h$ | 560.7 | 26.2 | 17.2 | 17.2
$\nu W^{\pm}Z$ | 492.3 | 12.3 | 9.3 | 9.3
Bkg. | - | - | - | 160.7

Table 3: Cut flow for the SM $\nu_{e}\bar{\nu}_{\mu}hh$ di-Higgs signal and
backgrounds cross sections at $e\mu$ collider with
$E_{e}=100\operatorname{GeV}$, $E_{\mu}=2.5\operatorname{TeV}$ (left) and
$E_{e}=170\operatorname{GeV}$, $E_{\mu}=6\operatorname{TeV}$ (right).

$\sigma$ [ab] | pre-cuts |  $m_{bb}$ cut | HHCUT
---|---|---|---
$\nu\nu hh(\to bbbb)$ | 12.2 | 6.8 | 5.5
$\nu\nu hZ$ | 49.8 | 9.3 | 5.2
$\nu\nu ZZ$ | 66.1 | 8.3 | 2.8
$\nu W^{\pm}h$ | 36.8 | 2.8 | 1.4
$\nu W^{\pm}Z$ | 40.2 | 2.2 | 0.9
Bkg. | - | - | 10.4

$\sigma$ [ab] | pre-cuts |  $m_{bb}$ cut | HHCUT
---|---|---|---
$\nu\nu hh(\to bbbb)$ | 69.9 | 37.6 | 31.5
$\nu\nu hZ$ | 205.5 | 33.1 | 19.3
$\nu\nu ZZ$ | 298.1 | 30.6 | 10.4
$\nu W^{\pm}h$ | 143.6 | 9.9 | 6.6
$\nu W^{\pm}Z$ | 139.6 | 6.5 | 2.8
Bkg. | - | - | 39.1

###  III.2 Asymmetric Colliders

The electron-muon ($e\mu$) collider, initially proposed in Refs. Hou (1996);
Choi  _et al._ (1998); Barger  _et al._ (1997), has drawn attentions recently
Lu  _et al._ (2021). As a hybrid option, it is often considered to benefit
from both the clean environment of electron colliders and the high beam energy
of muon colliders. In this work, we use two benchmark colliding energies to
illustrate the potential of the $e\mu$ collider for Higgs coupling
measurements:

  * •

A conventional 1 TeV setup with $E_{e}=100\operatorname{GeV}$ and
$E_{\mu}=3\operatorname{TeV}$,

  * •

A maximal colliding energy of 2 TeV with $E_{e}=170\operatorname{GeV}$ and
$E_{\mu}=6\operatorname{TeV}$.

The di-Higgs signal at the $e\mu$ collider is $4b+$$\not\mathrel{E}_{T}$,
similar to that at the $e^{+}e^{-}$ collider, and the background processes are
also comparable with the $e^{+}e^{-}$ collider and the
$\mu^{+}\mu^{-}$collider. However, due to non-zero lepton numbers in the
initial state, the $e\mu$ collider is free from the SM $s$-channel
backgrounds. Besides, there are no backgrounds from $ZZ/Zh$ production at the
$e\mu$ collider, and no need for a recoil mass cut. The collider simulation is
similar to that of the $e^{+}e^{-}$ collider. The calorimeter coverage range
of the detector is set as $-4<\eta<6$ to simulate the energy asymmetric
behavior of the $e\mu$ collider, where $\eta$ is the rapidity of the final
state objects. The following pre-selection cuts (pre-cuts) are applied:

$\displaystyle
n^{\ell}(p_{T}>10\operatorname{GeV})=0,~{}\not\mathrel{E}_{T}>10\operatorname{GeV},$
$\displaystyle p_{T}^{jet}>15\operatorname{GeV},-3.0<\eta^{jet}<5.0,~{}\Delta
R^{mn}>0.5.$ (10)

We then apply the same invariant mass cut and HHCUT as used for the symmetric
colliders, with the cut flows displayed in Table 3. With an integrated
luminosity of $3\operatorname{ab}^{-1}$, the signal strength of the total
cross section for di-Higgs production, $\mu_{hh}=\sigma_{hh}/\sigma_{hh}^{\rm
SM}$, can be constrained to $0.38<\mu_{hh}<1.78$ at the $1\operatorname{TeV}$
$e\mu$ collider and $0.75<\mu_{hh}<1.25$ at the $2\operatorname{TeV}$ $e\mu$
collider at 95% confidence level.

##  IV Analysis

Equipped with the simulation results, we explore the discovery potential of
lepton colliders on Higgs-gauge coupling, equivalently, $\kappa_{hhWW}$ and
$\kappa_{hWW}$. As stated before, to measure $\kappa_{hhWW}$ from the di-Higgs
production, the $\kappa_{hWW}$ must first be measured from the single Higgs
production. However, with a fixed value of $\kappa_{hWW}$, the total di-Higgs
cross section is a quadratic function of $\kappa_{hhWW}$; see Eq. (2).
Therefore, the measurement of $\kappa_{hhWW}$ may exhibit degeneracy. The
degeneracy can be resolved by analyzing the $m_{hh}$ distribution. As shown in
Fig. 5, the $m_{hh}$ distribution for different $\kappa_{hhWW}$ are
significantly different, which is a consequence of the unitarity violation
when $\kappa_{hhWW}\neq 1$; see Appendix A for details.

Figure 5: Normalized $m_{hh}$ distribution at 10 TeV $\mu^{+}\mu^{-}$collider
for different $\kappa_{hhWW}$, assuming $\kappa_{h}$ and $\kappa_{hWW}$ to be
the SM value.

Figure 6: Parameter space of $\kappa_{hWW}$ and $\kappa_{hhWW}$ that is
consistent with SM ($\kappa_{hWW}$ =$\kappa_{hhWW}$ =1) at 95% C.L. The left
panel gives the results from $1\operatorname{TeV}$ $e^{+}e^{-}$ collider and
$10\operatorname{TeV}$ $\mu^{+}\mu^{-}$collider, the right panel is the case
of $e\mu$ collider with $\sqrt{s}=1,2\operatorname{TeV}$. The dashed contour
in left panel stands for the discovery capability at the 10 TeV
$\mu^{+}\mu^{-}$collider.

In the following, we perform a binned analysis which utilizes the information
of $m_{hh}$ distribution. The likelihood function for a Poisson distribution
is given by

$L(\vec{\kappa})=\frac{\left(n_{\rm th}(\vec{\kappa})\right)^{n_{\rm
obs}}e^{-n_{\rm th}(\vec{\kappa})}}{n_{\rm obs}!},$ (11)

where $\vec{\kappa}=(\kappa_{hWW},\kappa_{hhWW})$ stands for the parameters in
the theory, $n_{\rm th}(\vec{\kappa})$ is the events number predicted by the
theory and $n_{\rm obs}$ is the observed events number. For the di-Higgs boson
production, the events number from theory prediction is $n_{\rm
th}(\vec{\kappa})=\left[\sigma_{hh}(\vec{\kappa})+\sigma_{\rm
bkg}\right]\mathcal{L}$ with an integrated luminosity of $\mathcal{L}$, where
$\sigma_{hh}$ is the di-Higgs production cross section after applying
selection cuts and $\sigma_{\rm bkg}$ is corresponding backgrounds. The
logarithm of likelihood function can be defined as

$\displaystyle-2\ln\lambda(\vec{\kappa}_{1}|\vec{\kappa}_{0})=-2\ln\frac{L(\vec{\kappa}_{1})}{L(\vec{\kappa}_{0})},$
(12)

which depicts the hypothesis $\vec{\kappa}_{1}$ is excluded versus
$\vec{\kappa}_{0}$ with
$\sqrt{-2\ln\lambda(\vec{\kappa}_{1}|\vec{\kappa}_{0})}\sigma$ confidence
level. For binned events, the likelihood function in Eq. (11) is modified as

$L(\vec{\kappa})=\prod_{i}^{N}\frac{\left(n^{i}_{\rm
th}(\vec{\kappa})\right)^{n^{i}_{\rm obs}}e^{-n^{i}_{\rm
th}(\vec{\kappa})}}{n^{i}_{\rm obs}!},$ (13)

where $n_{\rm th}^{i}(\vec{\kappa})$ and $n_{\rm obs}^{i}$ are the theory
predicted and observed event numbers in the $i^{\rm th}$ bin.

With the likelihood presented above, we discern the potential for measuring
Higgs-gauge-boson couplings. The region of
$-2\ln\lambda(\vec{\kappa}|\vec{\kappa}_{\rm SM})\leq 4$ is depicted in Fig.
6. Figure 6(a) displays results from the $e^{+}e^{-}$ collider and the
$\mu^{+}\mu^{-}$collider. The green region represents the
$\kappa_{hhWW}-\kappa_{hWW}$ parameter space that is consistent with the SM at
$95\%$ confidence level (C.L.) at the 1 TeV $e^{+}e^{-}$ collider, while the
orange region denotes the same for the 10 TeV $\mu^{+}\mu^{-}$collider. The
red star marks the SM prediction where $\kappa_{hWW}$ =$\kappa_{hhWW}$ =1. The
dashed contour in Fig. 6(a) stands for the $5\sigma$ discovery capability at
the 10 TeV $\mu^{+}\mu^{-}$collider. The blue band denotes the 95% C.L.
expected constraint on $\kappa_{hWW}$ from the future electron-proton
collider, namely $0.96<\kappa_{hWW}<1.04$ Sharma and Shivaji (2022), while the
current limit on $\kappa_{hWW}$ at the 13 TeV LHC is $0.93<\kappa_{hWW}<1.22$
Sirunyan  _et al._ (2019). Figure 6(b) presents results from asymmetric
colliders, specifically the $e\mu$ collider with $\sqrt{s}=1$ TeV (the green
region) and $\sqrt{s}=2$ TeV (the orange region).

From these figures, it is evident that the measurements become more accurate
with increasing collision energy due to the rise in event numbers with
$\sqrt{s}$, regardless of symmetric or asymmetric colliders. For example, in
the $e^{+}e^{-}$ collider and $e\mu$ collider with $\sqrt{s}=1$ TeV, the
signal and background event numbers are similar, as seen in Table 3.
Consequently, it is expected that the potentials are comparable, as shown by
the green range in Fig. 6(a) and (b).

As shown in Fig. 6, with the expected limit on $\kappa_{hWW}$, we obtained the
expected discovery sensitivity on $\kappa_{hhWW}$ as follows:

  * •

$5\sigma$ discovery region $\kappa_{hhWW}$ $<0.86$ or $\kappa_{hhWW}$ $>1.32$
at the 10 TeV $\mu^{+}\mu^{-}$collider,

  * •

$2\sigma$ exclusion region $\kappa_{hhWW}$ $<0.88$ or $\kappa_{hhWW}$ $>1.14$
at the 10 TeV $\mu^{+}\mu^{-}$collider,

  * •

$2\sigma$ exclusion region $\kappa_{hhWW}$ $<0.72$ or $\kappa_{hhWW}$ $>1.70$
at the 1 TeV $e^{+}e^{-}$ collider,

  * •

$2\sigma$ exclusion region $\kappa_{hhWW}$ $<0.72$ or $\kappa_{hhWW}$ $>1.66$
at the 1 TeV $e\mu$ collider,

  * •

$2\sigma$ exclusion region $\kappa_{hhWW}$ $<0.84$ or $\kappa_{hhWW}$ $>1.32$
at the 2 TeV $e\mu$ collider.

##  V Conclusion

The quartic Higgs-gauge-boson coupling $g_{hhVV}$ is essential to the
electroweak symmetry breaking mechanism and gauge boson mass generation, and
it can be measured in the VBF-type di-Higgs production. However, at hadron
colliders, the process is typically overshadowed by the gluon-fusion di-Higgs
production, making the extraction of $g_{hhVV}$ challenging. In contrast,
lepton colliders offer an ideal platform for probing the VBF-type di-Higgs
production.

In the work we explored the potential of lepton collider on the quartic Higgs-
gauge-boson coupling, or equivalently, its ratio to SM prediction
$\kappa_{hhWW}$. Detailed collider simulations including parton showering
effects were performed at the symmetric $e^{+}e^{-}$ collider,
$\mu^{+}\mu^{-}$collider and the asymmetric $e\mu$ collider. We found that the
collider’s potential to constrain $\kappa_{hhWW}$ primarily depends on the
collision energy. With an integrated luminosity of 3 ab-1, the $2\sigma$
exclusion region is $\kappa_{hhWW}$ $<0.72$ or $\kappa_{hhWW}$ $>1.70$ and
$\kappa_{hhWW}$ $<0.72$ or $\kappa_{hhWW}$ $>1.66$ at the 1 TeV $e^{+}e^{-}$
collider and $e\mu$ collider
($E_{e}=100\operatorname{GeV},E_{\mu}=3\operatorname{TeV}$), respectively. At
the 2 TeV $e\mu$ collider
($E_{e}=170\operatorname{GeV},E_{\mu}=6\operatorname{TeV}$), the $2\sigma$
exclusion region is $\kappa_{hhWW}$ $<0.84$ or $\kappa_{hhWW}$ $>1.32$, and it
is $\kappa_{hhWW}$ $<0.88$ or $\kappa_{hhWW}$ $>1.14$ at 10 TeV
$\mu^{+}\mu^{-}$collider. In addition we obtain the $5\sigma$ discovery region
of $\kappa_{hhWW}$ $<0.86$ or $\kappa_{hhWW}$ $>1.32$ at the 10 TeV
$\mu^{+}\mu^{-}$collider.

Note added: Near the completion of this work, Ref. Dávila  _et al._ (2023)
appeared with related content.

###### Acknowledgements.

The work of QHC, KC and XRW is partly supported by the National Science
Foundation of China under Grant No. 11725520 and No. 12235001. The work of YDL
is partly supported by the National Science Foundation of China under Grant
No. 12075257.

##  Appendix A Perturbative Unitarity

The nature of the VBF process can be ascertained by the factorized subprocess
of $W^{-}W^{+}\to hh$. This is executed through the Effective-$W$
approximation (EWA) method Dawson (1985) and the cross section is expressed as
the equation below:

$\sigma(\ell^{-}\ell^{+}\to\nu_{\ell}\bar{\nu}_{\ell}hh)=\sum_{i,j}\int
dx_{1}dx_{2}P^{i}_{W^{-}/\ell^{-}}(x_{1})P^{j}_{W^{+}/\ell^{+}}(x_{2})~{}\sigma(W^{-}_{i}W^{+}_{j}\to
hh).$ (14)

Here $P^{i}_{W^{-}/\ell^{-}}(x_{1})$ denotes the gauge boson $W^{-}$’s
probability distribution with polarization $i$ in the electron;
$x_{1}=(E_{\ell^{-}}-E_{\nu_{\ell}})/E_{\ell^{-}}$ is the energy fraction of
$\ell^{-}$ carried by $W^{-}$. The splitting functions for the collinear
radiation of real $W$-boson at leading order are Barger  _et al._ (1990).

$\displaystyle
P_{W/\ell}^{+}(x)=\frac{g^{2}}{32\pi^{2}x}\ln\left(\frac{4x^{2}E_{\ell}^{2}}{M_{W}^{2}}\right),$
$\displaystyle
P_{W/\ell}^{-}(x)=\frac{g^{2}}{32\pi^{2}x}(1-x)^{2}\ln\left(\frac{4x^{2}E_{\ell}^{2}}{M_{W}^{2}}\right),$
$\displaystyle P_{W/\ell}^{L}(x)=\frac{g^{2}}{16\pi^{2}x}(1-x).$ (15)

The matrix element of the subprocess $W^{\mu}(q_{1})W^{\nu}(q_{2})\to
h(k_{1})h(k_{2})$ is

$\displaystyle M^{\mu\nu}=$
$\displaystyle\left[\frac{3m_{h}^{2}g_{hWW}}{(\hat{s}-m_{h}^{2})v}+2g_{hhWW}+\left(\frac{g_{hWW}^{2}}{\hat{t}-m_{W}^{2}}+\frac{g_{hWW}^{2}}{\hat{u}-m_{W}^{2}}\right)\right]g^{\mu\nu}-$
$\displaystyle\frac{g_{hWW}^{2}}{m_{W}^{2}}\left[\frac{(q_{1}-k_{1})^{\mu}(q_{1}-k_{1})^{\nu}}{\hat{t}-m_{W}^{2}}+\frac{(q_{1}-k_{2})^{\mu}(q_{1}-k_{2})^{\nu}}{\hat{u}-m_{W}^{2}}\right],$
(16)

where $\hat{s}$ ($\hat{t}$, $\hat{u}$) are Mandelstam variables of the
subprocess, and $g_{hWW}$ and $g_{hhWW}$ are gauge couplings of the Higgs
boson with $W^{\pm}$ boson. The SM presents the couplings in the following:

$\displaystyle g_{hhWW}^{\rm SM}=\frac{m_{W}^{2}}{v^{2}},\quad g_{hWW}^{\rm
SM}=2\frac{m_{W}^{2}}{v},$ (17)

where $v=246\operatorname{GeV}$ refers to the electroweak symmetry breaking
scale. The couplings are parametrized as $g_{hWW}=\kappa_{hWW}g_{hWW}^{\rm
SM}$ and $g_{hhWW}=\kappa_{hhWW}g_{hhWW}^{\rm SM}$. The cross section of the
subprocess $W^{+}W^{-}\to hh$ can then be parametrized as follows:

$\displaystyle\sigma(W^{+}W^{-}\to hh)=$
$\displaystyle\kappa_{hhWW}^{2}\sigma_{a}+\kappa_{hWW}^{4}\sigma_{b}+\kappa_{hWW}^{2}\sigma_{c}+\kappa_{hhWW}\kappa_{hWW}^{2}\sigma_{ab}+$
$\displaystyle\kappa_{hhWW}\kappa_{hWW}\sigma_{ac}+\kappa_{hWW}^{3}\sigma_{bc},$
(18)

where $\sigma_{a}$, $\sigma_{b}$, and $\sigma_{c}$ parametrize the
contribution from Feynman diagrams (a), (b), and (c) depicted in Fig. 2,
whereas $\sigma_{ab}$, $\sigma_{ac}$, and $\sigma_{bc}$ are the corresponding
interferential contributions. The cross section of the complete process, i.e.,
$\ell^{-}\ell^{+}\to\nu_{\ell}\bar{\nu}_{\ell}hh$, can be parametrized
similarly. The dependence of the cross section of $W^{+}W^{-}\to hh$ on the
couplings $\kappa_{hWW}$ and $\kappa_{hhWW}$ can be inferred from the
parameters in Eq. (A). As shown in Fig. 8, the contributions from Feynman
diagrams (a) and (b) in Fig. 7 and their corresponding interference components
dominate the cross section.

Figure 7: Feynman diagrams representing $VV\to hh$ ($V=W^{\pm}/Z$) processes.
These correspond to the factorized subprocess shown in Fig. 2. Figure 8: The
energy dependence of each Feynman diagram contribution of the subprocess
$W^{+}_{L}W^{-}_{L}\to hh$. The solid line denotes the contribution from
Feynman diagram (a), (b) and (c), respectively; and the dashed line the
interference contribution.

In high-energy limits, the amplitude can be depicted by:

$\displaystyle\mathcal{M}_{LL}=M^{\mu\nu}\epsilon^{\mu}_{L}(q_{1})\epsilon^{\nu}_{L}(q_{2})=\frac{\hat{s}}{v^{2}}\big{(}\kappa_{hhWW}-\kappa_{hWW}^{2}\big{)}+\mathcal{O}(\hat{s}^{0}),$
(19)

for the longitudinally polarized vector bosons. This notion holds true as the
longitudinal polarization vector $\epsilon^{\mu}_{L}(q)$ approximates to
$q^{\mu}/m_{W}$. This expression reveals that high energy collision regions
primarily receive contributions from Feynman diagrams (a) and (b), a finding
that is consistent with the above numeric result. Furthermore, the violation
of the unitarity bound is inevitable unless the sum rule detailed below is
fulfilled:

$\displaystyle\kappa_{hhWW}-\kappa_{hWW}^{2}=0.$ (20)

Perturbative unitarity requires the partial wave amplitude $a_{l}$ Lee  _et
al._ (1977a, b) to satisfy

$\displaystyle|\operatorname{Re}{a_{l}}|\leq 1/2,$ (21)

with the definition of

$\displaystyle
a_{l}=\frac{1}{32\pi}\int_{-1}^{1}{\mathcal{M}}(\hat{s},\hat{t}){\text{P}}_{l}(\cos\theta)\operatorname{d}\cos\theta,$
(22)

and $\mathcal{M}(\hat{s},\hat{t})$ is amplitude of the process $W^{-}W^{+}\to
hh$ given in Eq. (19). Consequently, the unitarity bound constrains the
parameters $\kappa_{hWW}$ and $\kappa_{hhWW}$ for a given collision energy
$\sqrt{\hat{s}}$ of the subprocess; see Fig. 9. The colored regions in the
figure represent the unitarity-consistent areas for various collision
energies, while the solid black line indicates the sum rule in Eq. (20). The
figure demonstrates that, as $\sqrt{\hat{s}}=m_{hh}$ increases, deviations
from the sum rule must become less severe to maintain the theory consistence.

Figure 9: Parameter space of $\kappa_{hhWW}$ and $\kappa_{hWW}$ that satisfies
unitarity.

## References

  * Aad  _et al._ (2012) G. Aad _et al._ (ATLAS), Phys. Lett. B 716, 1 (2012), arXiv:1207.7214 [hep-ex] .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Phys. Lett. B 716, 30 (2012), arXiv:1207.7235 [hep-ex] .
  * Agashe  _et al._ (2005) K. Agashe, R. Contino, and A. Pomarol, Nucl. Phys. B 719, 165 (2005), arXiv:hep-ph/0412089 .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C 79, 421 (2019), arXiv:1809.10733 [hep-ex] .
  * Sharma and Shivaji (2022) P. Sharma and A. Shivaji, JHEP 10, 108 (2022), arXiv:2207.03862 [hep-ph] .
  * de Florian and Mazzitelli (2013) D. de Florian and J. Mazzitelli, Phys. Rev. Lett. 111, 201801 (2013), arXiv:1309.6594 [hep-ph] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), JHEP 07, 108 (2020), [Erratum: JHEP 01, 145 (2021), Erratum: JHEP 05, 207 (2021)], arXiv:2001.05178 [hep-ex] .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), Phys. Rev. D 108, 052003 (2023), arXiv:2301.03212 [hep-ex] .
  * Roloff  _et al._ (2020) P. Roloff, U. Schnoor, R. Simoniello, and B. Xu (CLICdp), Eur. Phys. J. C 80, 1010 (2020), arXiv:1901.05897 [hep-ex] .
  * Han  _et al._ (2021a) T. Han, D. Liu, I. Low, and X. Wang, Phys. Rev. D 103, 013002 (2021a), arXiv:2008.12204 [hep-ph] .
  * Han  _et al._ (2021b) T. Han, Y. Ma, and K. Xie, Phys. Rev. D 103, L031301 (2021b), arXiv:2007.14300 [hep-ph] .
  * Lu  _et al._ (2021) M. Lu, A. M. Levin, C. Li, A. Agapitos, Q. Li, F. Meng, S. Qian, J. Xiao, and T. Yang, Adv. High Energy Phys. 2021, 6693618 (2021), arXiv:2010.15144 [hep-ph] .
  * Accettura  _et al._ (2023) C. Accettura _et al._ , Eur. Phys. J. C 83, 864 (2023), arXiv:2303.08533 [physics.acc-ph] .
  * ILC (2013) (2013), arXiv:1306.6352 [hep-ph] .
  * Asner  _et al._ (2013) D. M. Asner _et al._ , in  _Community Summer Study 2013: Snowmass on the Mississippi_ (2013) arXiv:1310.0763 [hep-ph] .
  * Dong  _et al._ (2018) M. Dong _et al._ (CEPC Study Group), (2018), arXiv:1811.10545 [hep-ex] .
  * Aryshev  _et al._ (2022) A. Aryshev _et al._ (ILC International Development Team), (2022), arXiv:2203.07622 [physics.acc-ph] .
  * Budker (1969) G. I. Budker, Conf. Proc. C 690827, 33 (1969).
  * Parkhomchuk and Skrinsky (1983) V. V. Parkhomchuk and A. N. Skrinsky, Conf. Proc. C 830811, 485 (1983).
  * Neuffer (1983) D. Neuffer, Part. Accel. 14, 75 (1983).
  * Palmer  _et al._ (1996) R. Palmer _et al._ , Nucl. Phys. B Proc. Suppl. 51, 61 (1996), arXiv:acc-phys/9604001 .
  * Stratakis  _et al._ (2022) D. Stratakis _et al._ , (2022), arXiv:2203.08033 [physics.acc-ph] .
  * Abbott  _et al._ (2022) B. Abbott _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.08135 [hep-ex] .
  * Aimè  _et al._ (2022) C. Aimè _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07256 [hep-ph] .
  * Han  _et al._ (2022a) T. Han, Z. Liu, L.-T. Wang, and X. Wang, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07351 [hep-ph] .
  * Han  _et al._ (2022b) T. Han, T. Li, and X. Wang, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.05484 [hep-ph] .
  * Delahaye  _et al._ (2019) J. P. Delahaye, M. Diemoz, K. Long, B. Mansoulié, N. Pastrone, L. Rivkin, D. Schulte, A. Skrinsky, and A. Wulzer, (2019), arXiv:1901.06150 [physics.acc-ph] .
  * Long  _et al._ (2021) K. Long, D. Lucchesi, M. Palmer, N. Pastrone, D. Schulte, and V. Shiltsev, Nature Phys. 17, 289 (2021), arXiv:2007.15684 [physics.acc-ph] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Sjostrand  _et al._ (2008) T. Sjostrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * Cacciari  _et al._ (2012) M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C 72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * de Favereau  _et al._ (2014) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .
  * Hou (1996) G. W.-S. Hou, Nucl. Phys. B Proc. Suppl. 51, 40 (1996), arXiv:hep-ph/9605204 .
  * Choi  _et al._ (1998) S. Y. Choi, C. S. Kim, Y. J. Kwon, and S.-H. Lee, Phys. Rev. D 57, 7023 (1998), arXiv:hep-ph/9707483 .
  * Barger  _et al._ (1997) V. D. Barger, S. Pakvasa, and X. Tata, Phys. Lett. B 415, 200 (1997), arXiv:hep-ph/9709265 .
  * Dávila  _et al._ (2023) J. M. Dávila, D. Domenech, M. J. Herrero, and R. A. Morales, (2023), arXiv:2312.03877 [hep-ph] .
  * Dawson (1985) S. Dawson, Nucl. Phys. B 249, 42 (1985).
  * Barger  _et al._ (1990) V. D. Barger, K.-m. Cheung, T. Han, and R. J. N. Phillips, Phys. Rev. D 42, 3052 (1990).
  * Lee  _et al._ (1977a) B. W. Lee, C. Quigg, and H. B. Thacker, Phys. Rev. Lett. 38, 883 (1977a).
  * Lee  _et al._ (1977b) B. W. Lee, C. Quigg, and H. B. Thacker, Phys. Rev. D 16, 1519 (1977b).