[1]Valentin Karam

1]Department of Electrical Engineering and Computer Science, Massachusetts
Institute of Technology, 77 Massachusetts Ave., Cambridge, 02139,
Massachusetts, USA

# Parameter extraction for a superconducting thermal switch (hTron) SPICE
model

<EMAIL_ADDRESS>Owen Medeiros  Tareq El Dandachi  Matteo Castellani  Reed
Foster  Marco Colangelo  Karl Berggren  [

###### Abstract

Efficiently simulating large circuits is crucial for the broader use of
superconducting nanowire-based electronics. However, current simulation tools
for this technology are not adapted to the scaling of circuit size and
complexity. We focus on the multilayered heater-nanocryotron (hTron), a
promising superconducting nanowire-based switch used in applications such as
superconducting nanowire single-photon detector (SNSPD) readout. Previously,
the hTron was modeled using traditional finite-element methods (FEM), which
fall short in simulating systems at a larger scale. An empirical-based method
would be better adapted to this task, enhancing both simulation speed and
agreement with experimental data.

In this work, we perform switching current and activation delay measurements
on 17 hTron devices. We then develop a method for extracting physical fitting
parameters used to characterize the devices. We build a SPICE behavioral model
that reproduces the static and transient device behavior using these
parameters, and validate it by comparing its performance to the model
developed in a prior work, showing an improvement in simulation time by
several orders of magnitude. Our model provides circuit designers with a tool
to help understand the hTron’s behavior during all design stages, thus
promoting broader use of the hTron across various new areas of application.

###### keywords:

superconducting, nanowire, heater-Tron, hTron, SPICE, model, cryotron

##  1 Introduction

The field of superconducting electronics has demonstrated significant benefits
in areas such as quantum computing [1], imaging [2], neuromorphic computing
[3], and digital logic [4]. Josephson junction (JJ)-based circuits dominate
the field, thanks to their fast operation speeds and low power dissipation
[5].

However, JJs have low signals, low normal resistance, and are sensitive to
magnetic fields. In contrast, superconducting nanowires can operate in noisy
environments, offer high gains and exhibit high impedance and fan-out [6, 7,
8]. Therefore, they present a complementary device to JJs for certain
applications. In particular, the capability of nanowire-based circuits to
drive high loads makes them compatible with complementary metal-oxide-
semiconductor (CMOS) technologies [9, 10, 11].

Multiple nanowire-based devices have been introduced in the past, such as the
nanocryotron (nTron) [6] and the heater-nanocryotron (hTron) [12]. The nTron
is a three-terminal device that uses a constriction at the gate input to
suppress a superconducting channel. It has a maximum demonstrated clock
frequency of $615.4\text{\,}\mathrm{MHz}$ [13], but suffers from leakage
current between the gate and channel. The hTron device — in its multiplanar
version — uses the Joule heating from a normal heater to suppress a
superconducting channel deposited beneath it. An oxide layer separates the
heater layer from the superconducting layer, which isolates the layers
galvanically while coupling them thermally. The hTron is easier to fabricate
than the nTron due to the absence of constriction, and does not present the
issue of leakage currents [12].

These two nanowire-based devices allowed the creation of a logic family [14],
a memory cell [15], and are particularly used in Superconducting nanowire
single-photon detector (SNSPD) arrays, both for pulse amplification and
readout [16, 17, 18, 19]. The nTron can also be used to translate SFQ pulses
to CMOS [9]. When a higher impedance is needed, an nTron and hTron can be
combined in an amplification stage, e.g., to drive an LED and communicate
between superconducting neurons [20, 21]. If used in larger and more complex
circuits, superconducting nanowires thus have the potential to complement
Josephson junctions (JJs) and thus enable new superconducting electronics
applications.

However, the scaling of these nanowire-based circuits is currently limited to
a few devices per circuit [19, 18]. We can partly explain this deficiency by
the lack of simulation tools for superconducting nanowires, which is crucial
to their development. Baghdadi et al. developed a 3D electrothermal model for
the hTron using finite-element modeling (FEM) techniques, modeling heat
exchanges inside the device by solving differential heat equations [12]. This
model, while helpful at the device level during early development stages,
cannot easily simulate large-scale circuits, mainly because of the
differential equations stiffness and the difference of scales in the geometry
of these layered devices. A simpler, less accurate 0D model has also been
developed, which was implemented in SPICE [21]. However, this simpler model is
slow to solve, and suffers from convergence issues. Moreover, both models
suffer from a lack of agreement with experimental data. Indeed, the heat
equations require various physical parameters to be approximated from
literature or experiment [12]. This means that the model has to be tuned and
optimized each time the geometry is changed. An arbitrary heater and channel
widths cannot simply be plugged in to get a nice agreement with measurements,
because many of the material’s thermal parameters are geometry-dependent
(e.g., boundary resistance or diffusion coefficient).

These issues can be solved by developing a physics-informed behavioral model.
Instead of focusing on the device’s microscopic physics and heat exchanges,
behavioral models fit experimental data using a minimal set of parameters.
This approach allows physics-informed behavioral models to be more simple due
to their empirical basis, while also being robust thanks to fitting equations
arising from phenomenological physics.

More precisely, we can divide the hTron response to a given heater current
into a static and a transient response. The static response is defined by the
channel switching current, whereas the device activation delay — delay between
the input of a heater pulse and the channel switching — leads to the modeling
of the transient response. To model the hTron behavior, the critical current
and activation delay dependence on heater current have to be measured.

In this paper, we characterize 17 hTron devices from a single wafer and
explain their static behavior using only two physical parameters. We introduce
a systematic approach to extract these parameters from measurements of
critical current as a function of heater current. Furthermore, we demonstrate
the correlation of these parameters with the heater and channel widths. The
transient response of the device, which depends on the heat flow from the
heater to the channel through the oxide layer, is also modeled to fit
experiments, allowing us to accurately simulate the device close to the
maximum operating speed. Finally, we assess the simulation speed and accuracy
of our method, comparing it to the previous electrothermal model by applying
our parameter extraction method on measurement data from the previous hTron
study [12]. We achieve improved agreements with published experimental data
and successfully replicate published results. Moreover, our simulation time is
lower by several order of magnitude, making our approach compatible for use by
circuit designers. This novel behavioral model is a first step into a broader
use of superconducting nanowire-based devices in more complex circuits.

##  2 Methods

Figure 1: Device geometry and typical circuit
a) Device geometry, layers stackup, and typical circuit to perform I-V
measurements using the 4-points measurement, which is further detailed on the
b) schematics. c) I-V curve showing the influence of heater current on the
switching current.

In this section, we outline our methodological approach to constructing,
characterizing and modeling hTron devices. First, we delve into the
fabrication techniques. Following this, we explain the electrical measurement
procedures applied to the hTron devices, providing insights into their
functional characteristics. Lastly, we describe our modeling techniques,
including the specifics of our SPICE implementation.

###  2.1 Fabrication

Here we present the hTron fabrication process we used, which is comparable to
the one introduced by Baghdadi et al [12]. The Figure 1a depicts the device
layout together with the layer stackup.

The hTron fabrication process starts with the deposition of a uniform layer of
Niobium nitride (NbN) layer with thickness
$d_{\text{c}}=$20\text{\,}\mathrm{nm}$$ onto a silicon substrate using RF-
biased sputtering [22]. Subsequently, $5\text{\,}\mathrm{nm}$ Ti and
$50\text{\,}\mathrm{nm}$ Au were deposited with a liftoff process to pattern
the pads, that serve as wire-bonding connections for the hTron ports.
Following this, the channel is patterned into the underlying NbN layer with
electron beam lithography (EBL) followed by reactive ion etching. A
$100\text{\,}\mathrm{nm}$ silicon oxide layer is deposited using a PECVD
technique to isolate the heater and channel layers. Finally, the heaters are
patterned with another liftoff process using the same Ti/Au thickness ratio.

In total, 9 distinct geometry combinations were patterned on two
$10\times$10\text{\,}\mathrm{mm}$$ dies from a single wafer. The nanowire
configurations varied in heater width $w_{\mathrm{h}}$ and channel width
$w_{\mathrm{c}}$ of $100\text{\,}\mathrm{nm}$, $500\text{\,}\mathrm{nm}$, and
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and all had lengths of
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$. However, one device with
$w_{\mathrm{h}}=$100\text{\,}\mathrm{nm}$$ and
$w_{\mathrm{c}}=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ was damaged while
measuring, resulting in a total of 17 functional devices. As seen in Figure
1a, each terminal of the hTron is connected to 2 different pads, known as a
4-point measurement setup.

###  2.2 Measurement Setup

This section details the experimental setup used to measure the relation
between switching current $I_{\mathrm{SW}}$ as a function of the heater
current $i_{\mathrm{H}}$, $I_{\mathrm{SW}}(i_{\mathrm{H}})$, as shown in
Figure 1a and Figure 1b. This relation explains the steady-state behavior of
the hTron device. For our experiments, we used an immersion probe inside a
liquid helium dewar [23], setting a substrate temperature of
$T_{\mathrm{SUB}}=$4.2\text{\,}\mathrm{K}$$.

The switching current is determined by gradually sweeping the channel current
both positively and negatively using an arbitrary waveform generator (AWG)
until the device switches in the resistive state, resulting in a detectable
voltage. The AWG waveform frequency is set to $10\text{\,}\mathrm{kHz}$ so
that the period would be much larger than the hotspot growth time. We record
the voltage as close to the channel as possible using the 4-points setup. The
dependence of the switching current on the heater current is then obtained by
biasing the heater with a DC current $i_{\mathrm{H}}$ while performing the
channel sweep. A typical result is shown as a current-voltage (I-V) plot in
Figure 1c, where we can see the switching current being gradually lowered by
the increase of heater current. For high heater currents, the channel
temperature approaches $T_{\mathrm{C}}$, and thus the channel is fully
suppressed. To compensate for the intrinsic stochastic nature of the switching
characteristics of the devices, each datapoint is the result of an average of
100 trials [24].

The channel switching current can also be modulated by increasing the entire
substrate temperature. To characterize this effect, we used a commercial
temperature controller with a heater located inside the immersion probe
itself. By doing so, the whole sample is globally heated, allowing the
characterization of the direct dependence between the channel’s switching
current and temperature.

To conclude, the evolution of the switching current was characterized as a
function of the substrate temperature and heater current. Both measurements
highlight different components of the device behavior. Throughout this paper,
we use the switching current density
$J_{\mathrm{SW}}=I_{\mathrm{SW}}/(w_{\mathrm{c}}\cdot d_{\text{c}})$ instead
of the switching current, making our measurements independent on the channel
width.

###  2.3 Modeling

This subsection presents our approach to modeling the hTron device, focusing
on developing a behavioral model grounded in empirical data. We describe the
fitting functions and parameters, which are critical for accurately
replicating and predicting the hTron device’s behavior.

####  2.3.1 Physics-informed behavioral model

We developed a physics-informed behavioral model based on a set of equations
and fitting parameters to fit the $I_{\mathrm{SW}}(i_{\mathrm{H}})$
experimental data that we collected. In our approach, we treat the hTron as a
four-terminal black-box, focusing on modeling its behavior rather than its
intricate physical properties, which we consider unknown.

The fitting parameters can typically vary from one device to another, and
across two devices with same geometry. However, we found correlations between
the fitting parameters and the device’s widths ($w_{\mathrm{c}}$ and
$w_{\mathrm{h}}$), allowing us to predict the behavior of a device from its
geometry.

####  2.3.2 Fitting functions

We translated the flow of events governing the hTron static behavior by using
two analytical expressions, the first one estimating the channel temperature
from a heater current input, and the second one predicting the switching
current at this temperature.

First we estimate the channel temperature from the heater current:

$t_{\mathrm{Ch}}(i_{\mathrm{H}})=\left[(T_{\mathrm{C}}^{4}-T_{\mathrm{SUB}}^{4})\cdot\left(\frac{i_{\mathrm{H}}}{I_{\mathrm{H,SUPP}}}\right)^{\eta}+T_{\mathrm{SUB}}^{4}\right]^{\frac{1}{4}}$
(1)

$I_{\mathrm{H,SUPP}}$ is the suppressing current — the heater current at which
the channel is fully suppressed and reaches the critical temperature
$T_{\mathrm{C}}$. The parameter $\eta$ is the strength of the thermal
dependence between the heater current and the estimated channel temperature.
While we used $\eta=2$ in our model, its value will be discussed in Section 5.

We then find the unconstricted switching current density from the estimated
channel temperature, $\widehat{J}_{\mathrm{SW}}(t_{\mathrm{Ch}})$, defined as
the switching current density of the channel section located immediately below
the heater:

$\displaystyle\widehat{J}_{\mathrm{SW}}(t_{\mathrm{Ch}})=\widehat{J}_{\mathrm{C}}\cdot\left[1-\left(\frac{t_{\mathrm{Ch}}}{T_{\mathrm{C}}}\right)^{3}\right]^{2.1}\mathrm{,}\quad
t_{\mathrm{Ch}}\leq T_{\mathrm{C}}$ (2) $\displaystyle\text{with
}\widehat{J}_{\mathrm{C}}=\frac{\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})}{\left[1-\left(\frac{T_{\mathrm{SUB}}}{T_{\mathrm{C}}}\right)^{3}\right]^{2.1}}$

This function was introduced in [12]. The parameter $\widehat{J}_{\mathrm{C}}$
is the channel’s unconstricted critical current density, i.e., the switching
current at zero Kelvin of the channel section located immediately below the
heater. It can be obtained directly from
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$, the unconstricted switching
current density at substrate temperature.

By combining these two simple functions in the form
$\widehat{J}_{\mathrm{SW}}(t_{\mathrm{Ch}}(i_{\mathrm{H}}))=\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})$,
we were able to accurately fit our measurement data
$J_{\mathrm{SW}}(i_{\mathrm{H}})$. The fitting parameters
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ and $I_{\mathrm{H,SUPP}}$ are
extracted from measurements.

$\displaystyle\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})=\widehat{J}_{\mathrm{C}}\cdot\left[1-\left(\frac{1}{T_{\mathrm{C}}}\cdot\left[(T_{\mathrm{C}}^{4}-T_{\mathrm{SUB}}^{4})\cdot\left(\frac{i_{\mathrm{H}}}{I_{\mathrm{H,SUPP}}}\right)^{\eta}+T_{\mathrm{SUB}}^{4}\right]^{\frac{1}{4}}\right)^{3}\right]^{2.1}$
(3)

This equation models the switching current density of the channel portion
located directly below the heater. However, at low heater current the devices
often exhibit a constant switching current, or plateau. This plateau suggests
a defect or constriction along the channel in an area distant from the heater
(as detailed in Section 3: Results). In that case, the measured switching
current density at substrate temperature will differ from
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$. To match these observations,
the model switching current density can be rewritten as:

$\displaystyle
J_{\mathrm{MODEL}}(i_{\mathrm{H}})=\mathrm{min}\left\\{\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}}),J_{\mathrm{CONSTR}}\right\\},\;J_{\mathrm{CONSTR}}\leq\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$
(4)

The value $J_{\mathrm{CONSTR}}$ is the measured switching current when the
heater current is zero. In practice, it is equal to the measured switching
current density at substrate temperature, $J_{\mathrm{SW}}(T_{\mathrm{SUB}})$.
A value of $J_{\mathrm{CONSTR}}=\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$
would mean no plateau is observed. In practice, most of the measured devices
present a plateau at low heater currents.

Figure 2: End-to-end process for the hTron SPICE modeling.
a) Activation delay measurements allow to retrieve the transient device
response, and heater current-dependent switching current measurements define
the static device response. b) The heat transfer time constant
$\tau_{\mathrm{filter}}$ is computed from the activation delay. The fitting
parameters $I_{\mathrm{H,SUPP}}$ and
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$, together with the plateau level
$J_{\mathrm{CONSTR}}$ are extracted from $J_{\mathrm{SW}}(i_{\mathrm{H}})$
curves. c) The fitting parameters are further fed into SPICE, computing the
nanowire temperature, switching current, and retrapping current at each
simulation time step thanks to the according fitting functions. The heat
transfer through the oxide layer is simulated by a lumped-element RC circuit,
giving rise to a delay in the hTron response to an input heater pulse. d) The
switching and retrapping current are finally used to define the hotspot
behavior, that is either growing or retrapping depending on the current
flowing in the channel $I_{\mathrm{Ch}}$. The device geometry and NbN
resistance define the hotspot growth rate and normal resistance.

###  2.4 SPICE Implementation

In this subsection we describe the steps undertaken to build the simulation
model in SPICE [25], thus enabling the use of the modeled hTron device in
circuits. The SPICE model file is available in the supplementary material.

When a particular device with channel width $w_{\mathrm{c}}$ and heater width
$w_{\mathrm{h}}$ is to be characterized, the device geometry and wafer
material properties are gathered prior to performing electrical measurements.
The NbN thickness are measured using an ellipsometer [26], and the widths and
lengths are measured from scanning electron micrographs. The NbN sheet
resistance, which is used by the hotspot growth model, is also acquired. The
steps showing the building of a SPICE model from the electrical
characterization of a device are schematically represented in Figure 2.

In a first step (Figure 2a), electrical measurements must be performed. The
activation delay $\tau_{\mathrm{on}}$ (delay between the input of a heater
pulse and the channel switching) is measured to further define the transient
behavior of the device. The heater current-dependent switching current density
$J_{\mathrm{SW}}(i_{\mathrm{H}})$ is also measured to set the static device
behavior.

Subsequently (Figure 2b), the time constant of the heat transfer sub-circuit
$\tau_{\mathrm{filter}}$, the suppressing current $I_{\mathrm{H,SUPP}}$, the
unconstricted switching current at substrate temperature
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$, and the plateau level
$J_{\mathrm{CONSTR}}$ are extracted from the measurement data. Specifically,
$\tau_{\mathrm{filter}}$ is directly computed from $\tau_{\mathrm{on}}$, and
$I_{\mathrm{H,SUPP}}$, $\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$, and
$J_{\mathrm{CONSTR}}$ are obtained from $J_{\mathrm{SW}}(i_{\mathrm{H}})$
curves using the method detailed in Section 3 : Results.

In a third step (Figure 2c), the SPICE behavioral model uses the fitting
parameters to compute the switching current $I_{\mathrm{SW}}$ and retrapping
current $I_{\mathrm{R}}$ from a heater current input. The temporal behavior of
the heat transfer between the heater and the NbN layer through the oxide is
modeled using a simple heat transfer circuit with lumped elements (low-pass
filter), which time constant is $RC=\tau_{\mathrm{filter}}$. Without this
added sub-circuit, the temperature of the nanowire would start to grow
instantaneously with a heater current input, which does not reflect real
observations. The exact computation of $\tau_{\mathrm{filter}}$ can be found
in the supplementary material.

The channel temperature and switching current are then computed using
expression-based functions defined using the .func SPICE directive. This
directive computes the value of a custom expression without the need of adding
a voltage node, and a SPICE behavioral source allows one to probe the value of
the function. Due to the low computational cost and analytical nature of the
fitting functions, the SPICE implementation of our fitting functions is a
straight-forward process. This simplicity is an asset, as the solver updates
the temperature, switching current and retrapping current at each time step,
in response to the input heater current.

Finally, the channel current is compared to the switching and retrapping
currents to compute the hotspot resistance (Figure 2d). A channel current
above the switching current threshold would induce hotspot growth, after which
it would eventually reset by decreasing below the retrapping current. The
SPICE model of the non-linear inductor and hotspot behavior was based on the
superconducting nanowire SPICE model introduced by Berggren et al. [27],
however, the model was modified to embed a hotspot growth circuit based on the
built-in SPICE integrator function sdt(). The work by El Dandachi [28] has
proven this new approach to present multiple benefits over the previous
nanowire SPICE model.

##  3 Results

Figure 3: Measurement results and method for extracting model parameters
Experimental data of the switching current density of 9 different hTron
geometries, plotted against a) the substrate temperature and against b) the
heater current. Equation 2 was fitted to one of the data curves, showing a
smooth transition from $T_{\mathrm{SUB}}$ to $T_{\mathrm{C}}$. Heater-
dependent switching current curves show a plateau at low heater currents,
indicating that the heater current does not modulate the channel switching
current. c) Method developed to extract the model parameters from
$J_{\mathrm{SW}}(i_{\mathrm{H}})$ measurement plots. A line is fitted to the
linear part of the measurement data, which sets the non-physical fitting
parameters $\widetilde{I}_{\mathrm{H}}$ and $\widetilde{J}_{\mathrm{SW}}$
(indicated by $\Box$). The suppressing current $I_{\mathrm{H,SUPP}}$ and the
unconstricted switching current at substrate temperature
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ (indicated by $\triangle$), that
define the unconstricted model $\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})$,
are further derived from the non-physical parameters. In this example, a
constriction strongly limits the switching current of the measured device at
low heater current, but can be taken into account by defining
$J_{\mathrm{MODEL}}(i_{\mathrm{H}})=\mathrm{min}\\{\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}}),J_{\mathrm{CONSTR}}\\}$.

In this chapter we present the results of our research, going through the
fabrication, the devices characterization and the model design. We will
further analyze the model parameters and explain their dependence in device
geometry. Finally, we will compare the simulations and measurements of the
device activation delay.

###  3.1 Fabrication results

The deposited NbN film shows a critical temperature of
$T_{\mathrm{C}}\simeq$12.5\text{\,}\mathrm{K}$$, a thickness — measured with
ellipsometry — of $d_{\text{c}}=$23.6\text{\,}\mathrm{nm}$$ and a sheet
resistance of
$R_{\mathrm{SHEET}}=$77.9\text{\,}\mathrm{\SIUnitSymbolOhm}\mathrm{/}$\square$.
The length of the channel and heater strips is
$l_{\text{c}}=l_{\text{h}}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. Using
these parameters, a depairing current of
$J_{\text{dep}}=$16.44\text{\times}{10}^{10}\text{\,}\mathrm{A}\mathrm{/}\mathrm{m}\mathrm{{}^{2}}$$
was computed using the Equation 13 of Charaev et al [29]. Heater resistances
of $44.2\text{\,}\mathrm{\SIUnitSymbolOhm}$,
$10.4\text{\,}\mathrm{\SIUnitSymbolOhm}$, and
$7.2\text{\,}\mathrm{\SIUnitSymbolOhm}$ were measured on devices with heaters
widths of $0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ respectively, however, device-to-
device variations are likely to be observed. The heater resistance could be
engineered by tuning the heater length, which would not impair the device
operation. However, power is wasted to the substrate if the heater length is
greater than the channel width.

###  3.2 Device characterization

As shown on Figure 3a, the channel switching current density of measured
devices is smoothly reduced with temperature, until the channel is fully
suppressed at $T_{\mathrm{C}}\simeq$12.5\text{\,}\mathrm{K}$$. Moreover, the
experimental data agrees nicely with the fitting function
$J_{\mathrm{SW}}(T)=$4.48\text{\times}{10}^{10}$\cdot\left(1-\left(T/T_{\mathrm{C}}\right)^{3}\right)^{2.1}[\unit{/\squared}]$,
that was fitted to a particular measurement curve. All other curves are thus
expected to disagree with this fit.

On the other hand, Figure 3b shows the switching current density’s dependence
on heater current. This plot is used to extract the two fitting parameters
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ and $I_{\mathrm{H,SUPP}}$,
together with the additional value $J_{\mathrm{CONSTR}}$. Two comments can be
made on this plot : (1) a plateau appears for low heater currents, and (2) the
suppressing current is correlated to the heater width but not to the channel
width. These observations are discussed in Section 3.4 and Section 3.5,
respectively.

###  3.3 Parameter Extraction

Here we explain the developed step-by-step method to extract the two model
parameters from a particular $J_{\mathrm{SW}}(i_{\mathrm{H}})$ measurement
curve, represented schematically with a typical measurement curve on Figure
3c. The key idea is to recover the measurements that would have been observed
if the constriction outside the heated area (i.e., the plateau) would not have
been there. To do so, we used the information contained in the linear part of
the $J_{\mathrm{SW}}(i_{\mathrm{H}})$ measurement data:

  1. 1.

Fit a straight line to the linear part of the curve, ignoring the plateau on
the left side and the part on the right side (where the switching current
reaches zero).

  2. 2.

This linear fit defines $\widetilde{J}_{\mathrm{SW}}$ and
$\widetilde{I}_{\mathrm{H}}$, the intersection between the line and the y and
x axis, respectively.

  3. 3.

The fitting parameters can be recovered :
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})=\alpha\cdot\widetilde{J}_{\mathrm{SW}}$,
and $I_{\mathrm{H,SUPP}}=\beta\cdot\widetilde{I}_{\mathrm{H}}$. The $\alpha$
and $\beta$ constants are correction parameters obtained by fitting a straight
line to the fitting function directly. These constants simply account for the
fitting function’s curvature, and are valid for a same value of $\eta$. We
found $\alpha=0.88$ and $\beta=1.25$ for $\eta=2$. Plugging
$I_{\mathrm{H,SUPP}}$ into Equation 1 and
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ into Equation 2 gives the
unconstricted model $\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})$ as defined in
Equation 3, which represents the measurement we would observe if there was no
plateau.

  4. 4.

The final model expression — including the plateau — is defined in Equation 4:
$J_{\mathrm{MODEL}}(i_{\mathrm{H}})=\mathrm{min}\\{\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}}),J_{\mathrm{CONSTR}}\\}$,
the minimum between the unconstricted model
$\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})$ and the constriction level
$J_{\mathrm{CONSTR}}$. Any $J_{\mathrm{CONSTR}}$ value can be set in order to
set various constriction levels, depending on the application.

Figure 4: Results analysis: prediction of the model parameters
a) In black: measured switching current density at substrate temperature
$J_{\mathrm{SW}}(T_{\mathrm{SUB}})$ (which can be seen in the form of a
plateau) against the device channel width. In blue: constriction factor $C$,
defined as the measured critical current density divided by the predicted
switching current density. The average constriction factor is similar for
$w_{\mathrm{c}}=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ but decreases significantly for
$w_{\mathrm{c}}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. b) Predicted
fitting parameter $\widetilde{I}_{\mathrm{H}}$ dependence on the heater width,
with a power-law fit. For both plots, each point is an average over all device
with indicated heater or channel width.

###  3.4 Analysis of constrictions outside of the heated area

The presence of a plateau at low heater current in Figure 3b suggests that
most devices’ switching currents are limited by one or multiple weak spots
located along the channel, away from the heated area. At low heater currents,
these weak spots switch at lower channel currents than the channel portion
below the heater. We therefore observe a lower switching current than that of
the heated region. As the heater current increases, a point is reached where
the heated part of the channel starts to switch at a current lower than that
of the weak spot. This marks a significant change, as the plateau abruptly
collapses. Beyond this point, with further increases in heater current, the
switching occurs directly under the heater, where the switching current is at
its minimum. This process continues until the channel is completely
suppressed.

Constrictions are common candidates for weak spots in superconducting
nanowires [30]. In Figure 4a (in black), we plotted the switching current
density of the weak spot $J_{\mathrm{SW}}(T_{\mathrm{SUB}})$ against the
channel width. It can be seen that the switching current density of the weak
spot is similar for $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$-wide channels
and $0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$-wide channels, but decreases
significantly for $0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$-wide channels,
while also showing greater variance. The same trend can also be seen with the
constriction factor $C$, plotted in blue, defined as the ratio between the
measured switching current at substrate temperature
$J_{\mathrm{SW}}(T_{\mathrm{SUB}})$ and the ideal predicted critical current
at substrate temperature $\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$. The
fact that the devices with the narrowest channels are more constricted would
potentially result from the channel line-width roughness, which degrades as
the width decreases [31].

Finally, if our fabrication process was perfect,
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ should in theory be the same for
all devices from a same wafer, and should not depend on the channel width. In
practice, despite some variance, the average
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})$ value is indeed comparable
across devices with different channel widths. This would imply that the
channel’s average switching current density is similar among devices, despite
showing localised constrictions.

###  3.5 Analysis and prediction of the suppressing current

The suppressing current parameter $I_{\mathrm{H,SUPP}}$ has smaller values for
narrower heaters, and can be expressed as a function of the heater width
$w_{\mathrm{h}}$, as shown in Figure 4b. The function
$\widetilde{I}_{\mathrm{H}}(w_{\mathrm{h}})=16.9\cdot w_{\mathrm{h}}^{2/3}$,
resulting from a power-law fitting, predicts the suppressing current for
different heater widths. As reducing $w_{\mathrm{h}}$ increases both the
resistance and the dissipated power per unit area, the suppressing current
decreases faster for narrower heaters. Here, the fitting parameter, valid for
an entire wafer, embeds information about the oxide thickness and heater’s
materials properties.

Figure 5: Prediction of the hTron maximum operating speed
a) Simplified schematics of the measurement setup of the hTron activation
delay $\tau_{\mathrm{on}}$, defined as the channel response time following the
application of a heater current pulse of magnitude $I_{\mathrm{H}}$. The
heater current pulse duration is much larger than $\tau_{\mathrm{on}}$. We
plotted the measured and simulated activation delay against b) the channel
bias current (for a constant heater current), and c) the heater current
amplitude (for a constant channel current). The graphs highlight the model’s
predictive capability across a large span of input currents, despite a
deviation from actual measurements at lower current levels, due the simplicity
of our heat transfer model. The device heater width and channel width are
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$, respectively.

###  3.6 hTron maximum operating speed

The transient response of the hTron is defined by the heat transfer from the
heater to the channel through the oxide. It limits the operating speed of the
hTron devices, resulting in a non-zero device activation delay,
$\tau_{\mathrm{on}}$. This delay, measured as the time needed to observe a
channel switch after the input of a heater pulse, is a key aspect of the
device behavior.

In Figure 5a, we show the measurement setup used to measure the activation
delay $\tau_{\mathrm{on}}$ at a specific $(I_{\mathrm{H}},I_{\text{Ch}})$
operating point. The results of sweeping the channel current while maintaining
a fixed heater current are shown in Figure 5b. Figure 5c illustrates the
results when the heater current is varied with a fixed channel current. In
order to reproduce this data in SPICE, we first choose an operating point
$(I_{\mathrm{H}}=$1455.1\text{\,}\mathrm{\SIUnitSymbolMicro
A}$,I_{\text{Ch}}=$280.3\text{\,}\mathrm{\SIUnitSymbolMicro A}$)$, and extract
the activation delay from the experiments
($\tau_{\mathrm{on}}=$11.85\text{\,}\mathrm{ns}$$). We further obtain the RC
time constant of our heat transfer sub-circuit to
$RC=\tau_{\mathrm{filter}}=$10.0\text{\,}\mathrm{ns}$$ directly from the
activation delay value (details can be found in the supplementary material).
Sweeping the input currents in the same way in SPICE as we did in the
experiments, we observe that the model is reliable around the operating point
$(I_{\mathrm{H}}=$1455.1\text{\,}\mathrm{\SIUnitSymbolMicro
A}$,I_{\text{Ch}}=$280.3\text{\,}\mathrm{\SIUnitSymbolMicro A}$)$ and for high
currents.

This simple and efficient method allowed us to simulate the temporal behavior
of the hTron device, and can be used to determine the maximum operating speed
in real circuits. While the model diverges from the measurements at low
currents, the accuracy of this result could be improved with added degrees of
freedom in the heat transfer sub-circuit.

Figure 6: Reproduction of Baghdadi et al. [12] results from an
$I_{\mathrm{SW}}(i_{\mathrm{H}})$ measurement.
a) Circuit used to simulate the hTron behavior in SPICE, identical to the one
used in Figure 3 of [12]. The power dissipated in the heater increases the
channel’s temperature $t_{\mathrm{Ch}}$, decreasing the channel switching
current $I_{\mathrm{SW}}$. As the channel switches to the normal phase, the
bias current $I_{\text{Ch}}$ is transferred to a
$50\text{\,}\mathrm{\SIUnitSymbolOhm}$ load $R_{\text{Load}}$, allowing the
channel to reset back to the superconducting state. SPICE simulation result of
the b) effective channel temperature and of the c) electrical behavior of the
hTron device. The channel temperature and switching current were computed
using Equation 1 and Equation 2, respectively. The device and simulation
parameters were given by [12], while the model parameters
$I_{\mathrm{H,SUPP}}$ and $J_{\mathrm{SW}}(T_{\mathrm{SUB}})$ were extracted
from the $I_{\mathrm{SW}}(i_{\mathrm{H}})$ plots using our method presented in
this work.

##  4 Model performance and comparison with previous attempt

In this section, we applied our parameter extraction method to heater-
dependent switching current measurements from Figure 5 of Baghdadi et al.
[12], successfully replicating the device behavior as shown in Figure 6. The
measured device switching current at substrate temperature is approximately
$J_{\mathrm{SW}}(T_{\mathrm{SUB}})\simeq$120\text{\,}\mathrm{\SIUnitSymbolMicro
A}$$, and the $J_{\mathrm{SW}}(i_{\mathrm{H}})$ curve does not present a
plateau at low heater current.

Our model closely matches the measurements but slightly overestimates the
measured $J_{\mathrm{SW}}(T_{\mathrm{SUB}})$. This discrepancy might be due to
channel constriction. From their measurements, we extracted the parameters:
$I_{\mathrm{H,SUPP}}=$80.9\text{\,}\mathrm{\SIUnitSymbolMicro A}$$ and
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})=$122.9\text{\,}\mathrm{\SIUnitSymbolMicro
A}$$. Applying our method to the previous 3D electrothermal model results
curve from Figure 5 of [12] resulted in
$I_{\mathrm{H,SUPP}}=$79.5\text{\,}\mathrm{\SIUnitSymbolMicro A}$$ and
$\widehat{J}_{\mathrm{SW}}(T_{\mathrm{SUB}})=$116.1\text{\,}\mathrm{\SIUnitSymbolMicro
A}$$, against the previous model that predicted
$J_{\mathrm{SW}}(T_{\mathrm{SUB}})=$113.9\text{\,}\mathrm{\SIUnitSymbolMicro
A}$$. We computed a retrapping current value of
$I_{\mathrm{R}}=$20.3\text{\,}\mathrm{\SIUnitSymbolMicro A}$$ using the
relation given in [27], which is comparable to Baghdadi et al.’s measured
value of approximately
$I_{\mathrm{R}}\approx$21\text{\,}\mathrm{\SIUnitSymbolMicro A}$$.

In our simulation (Figure 6a), the hTron, shunted by a
$50\text{\,}\mathrm{\SIUnitSymbolOhm}$ load and biased at
$100\text{\,}\mathrm{\SIUnitSymbolMicro A}$, switches with a $1$-$\unit{}$
heater pulse, diverting the current into the resistor. The channel then
reverts to its superconducting state, allowing current return with a time
constant determined by the channel inductance and load resistance. Notably,
our SPICE simulation (Figure 6b and Figure 6c) doesn’t account for the
temperature rise from Joule heating when the channel is in the normal state.
However, this effect is taken into account by our hotspot growth rate equation
[27], which makes our model valid. The pulse and hotspot resistance, shown in
Figure 6c, align with Baghdadi et al.’s results (Figure 3c [12]), validating
our model.

In order to simulate the device transient behavior, we matched the device
activation time with the one observed in Figure 3c of [12], estimate as
$\tau_{\mathrm{on}}\approx$750\text{\,}\mathrm{p}\mathrm{s}$$. This gives an
RC time constant of $\tau_{\mathrm{filter}}=$1.6\text{\,}\mathrm{ns}$$ for the
lumped-element sub-circuit that simulates heat transfer.

All simulation parameters used in [12] are listed here:
$w_{\mathrm{c}}=$600\text{\,}\mathrm{nm}$$,
$w_{\mathrm{h}}=$500\text{\,}\mathrm{nm}$$, channel inductance
$L_{\text{Ch}}=$500\text{\,}\nH$$, $T_{\mathrm{SUB}}=$3\text{\,}\mathrm{K}$$,
$T_{\mathrm{C}}=$8.4\text{\,}\mathrm{K}$$, NbN thickness
$d_{\text{c}}=$20\text{\,}\mathrm{nm}$$,
$R_{\mathrm{SHEET}}=$470\text{\,}\mathrm{\SIUnitSymbolOhm}\mathrm{/}$\square$.
Note that the heater’s thickness and resistance, embedded in our fitting
parameters, are unused by our model.

Finally, we compared both models’ execution speeds by running the published
simplified 0D electrothermal model in SPICE, solving heat equations with
behavioral sources as implemented by Castellani [21], and comparing it to our
model in a $250\unit{\nano}$-long simulation of the circuit described in
Figure 6a. We have set the SPICE parameter $\text{{reltol}}=10^{-6}$ for both
models, as suggested by El Dandachi [28], and gathered the simulation results
in Table 1. In a first simulation run, we varied the maximum simulation
timestep $\Delta t_{\mathrm{max}}$ from $1\text{\,}\mathrm{ps}$ to
$10\text{\,}\mathrm{ps}$ over 11 steps, and our model took less than eleven
seconds to run, while the previous model took eight hours. In a second run, we
varied the maximum timestep from $1\text{\,}\mathrm{fs}$ to
$50\text{\,}\mathrm{ps}$ across 48 steps, and our model took less than four
hours to run, while the previous model had to be halted after two days because
it was stuck at a significantly slow simulation speed in the order of
$\approx$1\text{\,}\mathrm{f}\mathrm{s}\mathrm{/}\mathrm{s}$$. At this speed,
it would take almost 3 days to complete the $250\text{\,}\mathrm{ns}$ of
simulation before setting a higher value of maximum timestep.

To conclude, the enhanced simulation speed of our model can be mainly
explained by two factors. First, the previous model required solving a greater
number of nodes due to its approach of using behavioral sources to solve the
heat transfer equations of the entire system. Second, unlike our model, it did
not incorporate our enhanced hotspot integrator sub-circuit, potentially
leading to a higher number of timepoints [28]. While our simulations were
conducted on a personal computer, and performance can vary between different
machines, our results are significant. We successfully reduced typical
simulation times from hours to mere minutes or even seconds, making such
simulations feasible in practice on a personal computer for the first time,
opening the door for their widespread use.

|  Run 1: $\Delta t_{\text{max}}=$1\text{\,}\mathrm{ps}$$…$10\text{\,}\mathrm{ps}$ |  Run 2: $\Delta t_{\text{max}}=$1\text{\,}\mathrm{fs}$$…$50\text{\,}\mathrm{ps}$
---|---|---
# Time-points | Simulation Time | # Time-points | Simulation Time
Previous model | 813,628 |  $8.0$ hours | N/A | N/A
Our model | 2,045 |  $10.8$ seconds | 1,189,652 |  $3.3$ hours
Table 1: Comparison of simulation speed between the previous model designed by
Baghdadi et al. [12] and our model for the same SPICE circuit but different
time resolutions.

##  5 Discussion

This section discusses the challenges of predicting weak spots in
superconducting nanowires, the impact of hidden parameters in our model, and
the prediction of the device activation delay and operating speed,
highlighting both the strengths and limitations of our approach.

In the current state of research, accurately predicting the location and
number of weak spots — or areas with reduced critical current — along a
channel remains unfeasible. Weak spots located away from the heated area are
responsible for the plateau, which is not part of the behavior of the device,
but rather due to channel width-dependent effects arising from the fabrication
process. Therefore, we choose to get rid of the plateau during the parameter
extraction, effectively fitting to the linear part of the
$J_{\mathrm{SW}}(i_{\mathrm{H}})$ plots. As a result, the ideal model
$\widehat{J}_{\mathrm{SW}}(i_{\mathrm{H}})$ is a hypothetical curve that
probes the critical current directly below the heater, theoretically
representing the channel’s critical current as if there was no weak spots. The
plateau can then be set to any constriction level $J_{\mathrm{CONSTR}}$. While
we observed that the plateau level is smaller for narrower channels, it cannot
be accurately predicted due to its intrinsic stochastic nature. Locating and
removing constrictions in superconducting nanowires would be beneficial both
in photon detection and device operation.

Regarding the parameter $\eta$ in Equation 1, we used the value of $\eta=2$ in
this paper based on its optimal fit with Baghdadi et al.’s measurement data
[12]. However, Butters suggests $\eta=3$ (See Equation 3.30 of [23]), implying
a stronger dependence of the channel temperature with the heater current. The
value of $\eta=3$ aligns more closely with our measurement data, especially
when $i_{\mathrm{H}}$ is close to $I_{\mathrm{H,SUPP}}$. However, our devices
are severely constricted, showing a plateau in most cases. Thus more
experimental data on less constricted devices would help determine the best
value for $\eta$. While a more complicated equation could potentially better
describe the channel temperature as a function of heater current, the strength
of our approach lies in its simplicity: introducing an additional free
parameter would counteract this benefit, making our model less straightforward
without a substantial increase in predictive accuracy for our specific
application.

In this study, the measured devices are from a single wafer, implying that
some geometrical parameters such as the oxide thickness do not explicitly
appear in the fitting equations, and are thus hidden. Consequently, new
measurements and recalibration of the model’s parameters would be required to
predict the behavior of devices with a different oxide thickness. Even though
the previous model does not present hidden parameters, its prediction
capabilities are limited as well. Indeed, it contains various thermal
parameters that are either extracted from literature or fitted from
measurements, and cannot be predicted from geometry. In contrast, our high-
level curve-fitting approach simplifies the analysis, allows for a systematic
process of parameters extraction for any geometry, and outperforms the
previous attempt in accuracy and speed.

Finally, the heat transfer from the heater to the channel through the oxide is
a critical factor limiting the operating speed of the devices, resulting in a
non-zero device activation delay, $\tau_{\mathrm{on}}$. In complicated
circuits, this delay becomes an increasingly significant figure for circuit
designers. Indeed, accurately predicting and simulating the maximum operating
speed of a real-world circuit is essential. In our approach, we intentionally
over-simplified the way we modeled the heat transfer, offering the advantage
of only including one parameter — an RC time constant denoted
$\tau_{\mathrm{filter}}$, which can be directly obtained from
$\tau_{\mathrm{on}}$ —, making it straight-forward to implement in practice.
However, as seen in Figure 5, our model’s prediction of the activation delay
is limited for low input currents. At higher currents, though, our model
accurately predicts the activation delay, allowing one to simulate the
temporal behavior of a hTron-based circuit close to the maximum operating
speed. As a final note, while the minimum measured activation delay is in the
order of $\qtyrange{1}{10}{}$ for our particular device, it is noteworthy that
the fabricated devices were not designed for speed. A thinner oxide layer, as
well as narrower heater wires would greatly improve the reaction speed of
hTron devices.

##  6 Conclusion

Our approach to modeling the hTron device relies on the fitting of simple
physics-informed equations to a relatively large amount of experimental data.
All the fitting parameters are extracted from basic measurements, and no
arbitrary parameter has to be tuned to better fit our measurements. The static
behavior was modeled from the critical current measurements performed on 17
hTrons having 9 different geometries. We were able to extract the two relevant
fitting parameters, and extrapolate them for hTron geometries that were not
measured. Moreover, we were able to simulate the hTron transient response over
a large span of current inputs thanks to activation delay measurements — delay
between the application of a heater current pulse and the switching of the
channel. This transient behavior is critical in circuit design as it sets the
maximum operating speed of a circuit. We applied our model on measurement data
published by Baghdadi et al. [12], and obtained a better agreement than their
attempt. Moreover, we compared the simulation speed of both approaches on a
personal computer in SPICE, and obtained a result in seconds, while the
previous approach completed the same simulation in hours. Consequently, our
model and SPICE implementation is tailored for efficient design iterations in
nanowire-based electronics. Finally, our model effectively simulates devices
for any constriction — or plateau — level, a critical feature given the
unpredictability of these occurrences during circuit design.

Our findings corroborate the positioning of the hTron as an alternative to the
nTron device for applications requiring both high output impedance and
electrical isolation between gate and channel. In a broader view point,
superconducting nanowires could complement Josephson Junctions in areas where
they are currently lacking. Looking forward, we will use our model to solve
larger and more complex circuit like SNSPD arrays readouts. The high
simulation speed capability, accuracy, and simplicity of our approach
positions this work as a promising first step for the building of useful
simulation tools for nanowire-based circuits.

## Acknowledgements

The initial stages of the research were sponsored by the U.S. Department of
Energy, Office of Science, Office of Basic Energy Sciences, under Award Number
DE-AC02-07CH11359. The completion of data analysis and presentation were
sponsored by the National Science Foundation under Grant No. OMA-2137723. The
authors thank Alessandro Restelli, Joshua Bienfang and Ilya Charaev for
helpful scientific discussions. V.K. would like to thank Edoardo Charbon from
the Advanced Quantum Architecture (AQUA) Laboratory at the Swiss Federal
Institute of Technology (EPFL). O.M. acknowledges support from the NDSEG
Fellowship program. M. Colangelo acknowledges support from MIT Claude E.
Shannon award.

## References

  *   * Ladd et al. [2010] Ladd, T.D., Jelezko, F., Laflamme, R., Nakamura, Y., Monroe, C., O’Brien, J.L.: Quantum computers. nature 464(7285), 45–53 (2010)
  * Oripov et al. [2023] Oripov, B.G., Rampini, D.S., Allmaras, J., Shaw, M.D., Nam, S.W., Korzh, B., McCaughan, A.N.: A superconducting-nanowire single-photon camera with 400,000 pixels. arXiv preprint arXiv:2306.09473 (2023)
  * Islam et al. [2023] Islam, M.M., Alam, S., Hossain, M.S., Roy, K., Aziz, A.: A review of cryogenic neuromorphic hardware. Journal of Applied Physics 133(7) (2023)
  * Likharev [1993] Likharev, K.K.: In: Weinstock, H., Ralston, R.W. (eds.) Rapid Single-Flux-Quantum Logic, pp. 423–452. Springer, Dordrecht (1993). https://doi.org/10.1007/978-94-011-1918-4_14 . https://doi.org/10.1007/978-94-011-1918-4_14
  * Soloviev et al. [2017] Soloviev, I.I., Klenov, N.V., Bakurskiy, S.V., Kupriyanov, M.Y., Gudkov, A.L., Sidorenko, A.S.: Beyond moore’s technologies: operation principles of a superconductor alternative. Beilstein journal of nanotechnology 8(1), 2689–2710 (2017)
  * McCaughan and Berggren [2014] McCaughan, A.N., Berggren, K.K.: A superconducting-nanowire three-terminal electrothermal device. Nano letters 14(10), 5748–5753 (2014)
  * Zheng et al. [2020] Zheng, K., Zhao, Q.-Y., Lu, H.-Y.-B., Kong, L.-D., Chen, S., Hao, H., Wang, H., Pan, D.-F., Tu, X.-C., Zhang, L.-B., et al.: A superconducting binary encoder with multigate nanowire cryotrons. Nano Letters 20(5), 3553–3559 (2020)
  * Huang et al. [2023] Huang, Y.-H., Zhao, Q.-Y., Chen, S., Hao, H., Wang, H., Guo, J.-W., Tu, X.-C., Zhang, L.-B., Jia, X.-Q., Chen, J., et al.: Splitter trees of superconducting nanowire cryotrons for large fan-out. Applied Physics Letters 122(9) (2023)
  * Zhao et al. [2017] Zhao, Q.-Y., McCaughan, A.N., Dane, A.E., Berggren, K.K., Ortlepp, T.: A nanocryotron comparator can connect single-flux-quantum circuits to conventional electronics. Superconductor Science and Technology 30(4), 044002 (2017)
  * Toomey et al. [2019] Toomey, E., Onen, M., Colangelo, M., Butters, B., McCaughan, A., Berggren, K.: Bridging the gap between nanowires and josephson junctions: a superconducting device based on controlled fluxon transfer. Physical review applied 11(3), 034006 (2019)
  * McCaughan et al. [2019] McCaughan, A.N., Verma, V.B., Buckley, S.M., Allmaras, J., Kozorezov, A., Tait, A., Nam, S., Shainline, J.: A superconducting thermal switch with ultrahigh impedance for interfacing superconductors to semiconductors. Nature electronics 2(10), 451–456 (2019)
  * Baghdadi et al. [2020] Baghdadi, R., Allmaras, J.P., Butters, B.A., Dane, A.E., Iqbal, S., McCaughan, A.N., Toomey, E.A., Zhao, Q.-Y., Kozorezov, A.G., Berggren, K.K.: Multilayered heater nanocryotron: A superconducting-nanowire-based thermal switch. Physical Review Applied 14(5), 054011 (2020)
  * Zheng et al. [2019] Zheng, K., Zhao, Q.-Y., Kong, L.-D., Chen, S., Lu, H.-Y.-B., Tu, X.-C., Zhang, L.-B., Jia, X.-Q., Chen, J., Kang, L., et al.: Characterize the switching performance of a superconducting nanowire cryotron for reading superconducting nanowire single photon detectors. Scientific reports 9(1), 16345 (2019)
  * Buzzi et al. [2023] Buzzi, A., Castellani, M., Foster, R.A., Medeiros, O., Colangelo, M., Berggren, K.K.: A nanocryotron memory and logic family. Applied Physics Letters 122(14) (2023)
  * Butters et al. [2021] Butters, B.A., Baghdadi, R., Onen, M., Toomey, E.A., Medeiros, O., Berggren, K.K.: A scalable superconducting nanowire memory cell and preliminary array test. Superconductor Science and Technology 34(3), 035003 (2021)
  * Oripov et al. [2023] Oripov, B., Rampini, D., Allmaras, J., Shaw, M., Nam, S.W., Korzh, B., McCaughan, A.: Thermally coupled imager with 400,000 pixels (2023)
  * Zheng et al. [2020] Zheng, K., Zhao, Q.-Y., Lu, H.-Y.-B., Kong, L.-D., Chen, S., Hao, H., Wang, H., Pan, D.-F., Tu, X.-C., Zhang, L.-B., et al.: A superconducting binary encoder with multigate nanowire cryotrons. Nano Letters 20(5), 3553–3559 (2020)
  * Foster et al. [2023] Foster, R.A., Castellani, M., Buzzi, A., Medeiros, O., Colangelo, M., Berggren, K.K.: A superconducting nanowire binary shift register. Applied Physics Letters 122(15) (2023)
  * Castellani et al. [2023] Castellani, M., Medeiros, O., Foster, R.A., Buzzi, A., Colangelo, M., Bienfang, J.C., Restelli, A., Berggren, K.K.: A nanocryotron ripple counter integrated with a superconducting nanowire single-photon detector for megapixel arrays. arXiv preprint arXiv:2304.11700 (2023)
  * Shainline et al. [2019] Shainline, J.M., Buckley, S.M., McCaughan, A.N., Chiles, J.T., Jafari Salim, A., Castellanos-Beltran, M., Donnelly, C.A., Schneider, M.L., Mirin, R.P., Nam, S.W.: Superconducting optoelectronic loop neurons. Journal of Applied Physics 126(4) (2019)
  * Castellani [2020] Castellani, M.: Design of superconducting nanowire-based neurons and synapses for power-efficient spiking neural networks. Master’s thesis, Politecnico di Torino (2020)
  * Dane et al. [2017] Dane, A.E., McCaughan, A.N., Zhu, D., Zhao, Q., Kim, C.-S., Calandri, N., Agarwal, A., Bellei, F., Berggren, K.K.: Bias sputtered nbn and superconducting nanowire devices. Applied Physics Letters 111(12) (2017)
  * Butters [2022] Butters, B.A.: Digital and microwave superconducting electronics and experimental apparatus. PhD thesis, Massachusetts Institute of Technology (2022)
  * Lin and Roy [2013] Lin, S.-Z., Roy, D.: Role of kinetic inductance in transport properties of shunted superconducting nanowires. Journal of Physics: Condensed Matter 25(32), 325701 (2013)
  * Kundert [2003] Kundert, K.: The Designer’s Guide to Spice and Spectre. Kluwer Academic Publishers, Boston (2003)
  * Medeiros et al. [2019] Medeiros, O., Colangelo, M., Charaev, I., Berggren, K.K.: Measuring thickness in thin nbn films for superconducting devices. Journal of Vacuum Science & Technology A 37(4) (2019)
  * Berggren et al. [2018] Berggren, K.K., Zhao, Q.-Y., Abebe, N., Chen, M., Ravindran, P., McCaughan, A., Bardin, J.C.: A superconducting nanowire can be modeled by using spice. Superconductor Science and Technology 31(5), 055010 (2018)
  * El Dandachi [2023] El Dandachi, T.: Efficient simulation of large-scale superconducting nanowire circuits. Master’s thesis, Massachusetts Institute of Technology (2023)
  * Charaev et al. [2017] Charaev, I., Silbernagel, T., Bachowsky, B., Kuzmin, A., Doerner, S., Ilin, K., Semenov, A., Roditchev, D., Vodolazov, D.Y., Siegel, M.: Proximity effect model of ultranarrow nbn strips. Physical Review B 96(18), 184517 (2017)
  * Kerman et al. [2007] Kerman, A.J., Dauler, E.A., Yang, J.K., Rosfjord, K.M., Anant, V., Berggren, K.K., Gol’tsman, G.N., Voronov, B.M.: Constriction-limited detection efficiency of superconducting nanowire single-photon detectors. Applied Physics Letters 90(10) (2007)
  * Colangelo et al. [2022] Colangelo, M., Walter, A.B., Korzh, B.A., Schmidt, E., Bumble, B., Lita, A.E., Beyer, A.D., Allmaras, J.P., Briggs, R.M., Kozorezov, A.G., et al.: Large-area superconducting nanowire single-photon detectors for operation at wavelengths up to 7.4 $\mu$m. Nano Letters 22(14), 5667–5673 (2022)


# Vacuum Stability vs. Positivity in Real Singlet Scalar Extension of the
Standard Model

Parsa Ghorbani  Dipartimento di Fisica dell’Università di Pisa, Italy
INFN, Sezione di Pisa, Italy

###### Abstract

We assume a generic real singlet scalar extension of the Standard Model living
in the vacuum $(v,w)$ at the electroweak scale with $v=246$ GeV and $w$ being
respectively the Higgs and the singlet scalar vacuum expectation values. By
requiring absolute vacuum stability for the vacuum $(v,w)$, the positivity
condition and the perturbativity up to the Planck scale, we show that the
viable space of parameters in the model is strongly constrained for various
singlet scalar vacuum expectation values $w=0.1,1,10,100$ TeV. Also, it turns
out that the singlet scalar mass can be from a few GeV up to less than TeV.

##  1 Introduction

The stability of the vacuum in the Standard Model (SM) was first utilized as
an implement to put a bound on the Higgs mass and the mass of fermions in the
SM framework[5]. After the discovery of the Higgs particle by ATLAS and CMS
experiments at the LHC in 2012, the value of the Higgs mass is determined
accurately to be around $125$ GeV [1, 6]. The mass of the top quark (as the
heaviest quark) was already known to be about $176$ GeV [2]. Having the Higgs
and the top quark masses in hand, and knowing the value of the Higgs vacuum
expectation value (VEV) to be $v=246$ GeV, the status of the vacuum stability
in the SM becomes lucid; the SM vacuum starts to be metastable at energy scale
around $10^{10}$ GeV [3, 8, 4]. This happens because in the SM, the top quark
has a large negative contribution in the renormalization group equations (RGE)
for the Higgs quartic coupling $\lambda_{\text{h}}$, so that the Higgs quartic
coupling becomes negative at higher energy scales which makes the vacuum with
$v\neq 0$ metastable.

There is a consensus in the literature that in the presence of more scalars,
the vacuum can become stable up to the Planck scale. For instance, the vacuum
stability in extensions of the SM by adding an extra real scalar with
$\mathbb{Z}_{2}$ symmetry is studied in [17, 11], employing a complex scalar
to address the vacumm stability in [9, 18], using scalars in scale invariant
extension of the SM in [13, 20], investigation of the di-Higgs production in
singlet scalar model in [7], studying the vacuum stability of 2HDM at the
electroweak (EW) scale in [12], and stablizing the vacuum by Higgs–inflaton
mixing in [10].

The point that we want to emphasize in this paper is the importance of the
positivity condition (i.e. the requirement of having a positive definite
potential in all scales), when studying the vacuum stability in a given model.
At the EW scale, say at the scale $\mathcal{O}(m_{t})$, the free parameters of
an extended SM model must be chosen in a way to respect the positivity
condition. However, when solving the RGEs there is no guarantee that the
positivity condition will be satisfied in higher energy scales. This should be
considered alongside the possible change of the vacuum structure at higher
scales.

In the case of the SM, the vacuum stability and the positivity are delicately
related. The potential in the SM is given by
$V_{\text{SM}}=-\mu_{\text{h}}h^{2}/2+\lambda_{\text{h}}h^{4}/4$ for which if
$\lambda_{\text{h}}>0$ the theory develops a non-zero VEV for the Higgs
scalar. If $\lambda_{\text{h}}<0$ the Higgs VEV can only be vanishing. Due to
the large negative contribution of the heavy top quark in the RGE for
$\lambda_{\text{h}}$, the Higgs quartic coupling becomes negative at the scale
$10^{10}$ and thereafter the Higgs non-zero VEV is no longer a minimum of the
theory. Therefore, in the SM the sign of the quartic coupling
$\lambda_{\text{h}}$ changes the structure of the vacuum and deals with the
vacuum stability. At the same time, the sign of the Higgs quartic coupling
confirms or violates the positivity of the potential. If
$\lambda_{\text{h}}<0$ the theory is no longer well defined before thinking
about the vacuum structure of the model. Therefore, in the case of the SM, the
quartic coupling $\lambda_{\text{h}}$ plays a dual role as a tuner for the
vacuum stability and the positivity of the model.

When more scalars are added in the theory, the vacuum stability and the
positivity condition should be investigated separately. Both conditions even
if consistent at the EW scale, may become in conflict at higher scales. This
might resemble although maybe not directly related to the situation that for a
random set of parameters in a multi-scalar potential having at the same time a
small Higgs mass and a small cosmological constant is not very probable (see
[16]). In general with more scalars the vacuum structure of the model gets
more complicated; as the number of vacua grows rapidly with the number of
scalars (see [14] for two-scalar example). The positivity condition which is
obtained from the quartic part of the potential, can be very involved in
general with more scalars [19], unless some symmetry is applied.

The perturbativity of the theory is another constraint which must be taken
into account when running the parameters of the model to higher scales. The
main question is how consistent the vacuum stability constraint, the
positivity condition and the perturbativity are in a multi-scalar theory from
the electroweak up to the Planck scale.

In this article, as the first step to answer the question posed above, we
extend the SM with a generic real scalar potential including also the scalar
cubic term and linear scalar-Higgs interaction. We will study the scale
evolution of the only vacuum of the model i.e. $(v,w)$ at the EW scale, and
will argue the absolute vacuum stability up to the Planck scale. By absolute
vacuum stability we mean that the potential possess only one single minimum
for whole range of the energy scale here from the EW up to the Planck. This
will be confronted with the scale evolution of the positivity condition as
well as the perturbativity up to the Planck scale. Using the aformentioned
constraints we put strong bound on free parameters of the model. For
simplicity we use the term SPP conditions when we consider the stability,
positivity and perturbativity conditions altogether.

The rest of the paper is arranged as the following. In section 2 we introduce
the model giving the details of the vacuum solution $(v,w)$ and the positivity
condition. In section 3 we will discuss the RGEs, and will present our
numerical analysis in section 4. We will bring a summery of the results in
section 5.

##  2 Vacua in Singlet Scalar Model

The vacuum stability in the real singlet scalar extension of the SM with
$\mathbb{Z}_{2}$ symmetry has been studied in [11]. The presence of the
$\mathbb{Z}_{2}$ symmetry simplifies the model considerably. For such model
provided that the positivity condition is taken into account, if the vacuum
solution $(v,w)$ is a minimum at a given scale $\Lambda$, it remains the
global minimum in all scales because $(v,w)$ and other extremum solutions of
the $\mathbb{Z}_{2}$ symmetric model i.e. $(0,0)$, $(v,0)$ and $(0,w)$ cannot
be minimum at the same time (see [15] on vacuum structure of the
$\mathbb{Z}_{2}$ symmetric singlet scalar model).

Here we consider instead a generic real scalar extension of the SM without the
$\mathbb{Z}_{2}$ symmetry,

$w$ |  $100$ GeV |  $1$ TeV |  $10$ TeV |  $100$ TeV
---|---|---|---|---
$\lambda_{\text{h}}$ | $(0.39,0.51)$ | $(0.1,0.54)$ | $(1.9\times 10^{-4},0.39)$ | $(0.95,1)$
$\lambda_{\text{s}}$ | $(0.25,1)$ | $(0.015,0.032)$ | $(2.4,4.07)\times 10^{-4}$ | $(1.03,1.13)\times 10^{-5}$
$\lambda_{\text{hs}}$ | $(-0.17,-0.02)$ | $(-0.033,0.015)$ | $(-2900,-4.53)\times 10^{-6}$ | $(-2.37,-0.96)\times 10^{-6}$
$\kappa_{\text{s}}$ | $(-1,1)$ | $(-1,1)$ | $(-1,1)$ | $(0.49,0.64)$
$\kappa_{\text{hs}}$ | $(-1,1)$ | $(-0.054,1)$ | $(-1,0.097)$ | $(-0.55,-0.098)$
$m_{s}$ |  $191$-$257$ GeV |  $218$-$308$ GeV |  GeV $185$-$334$ GeV |  $553$-$574$ GeV
Table 1: The allowed region for the couplings
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}}$,
$\kappa_{\text{s}},\kappa_{\text{hs}}$ and $m_{s}$ for different singlet
scalar VEV benchmarks $w=0.1,1,10,100$ TeV, respecting the absolute vacuum
stability and positivity condition for the vacuum $(v,w)$ at the electroweak
scale with the assumptions $m_{H}=125$ GeV, $v_{H}=246$ GeV and $m_{s}>m_{H}$.

$V(h,s)=-\frac{1}{2}\mu^{2}_{\text{h}}h^{2}+\frac{1}{4}\lambda_{\text{h}}h^{4}-\frac{1}{2}\mu^{2}_{\text{s}}s^{2}+\frac{1}{3}\kappa_{\text{s}}s^{3}+\frac{1}{4}\lambda_{\text{s}}s^{4}+\frac{1}{2}\kappa_{\text{hs}}h^{2}s+\frac{1}{4}\lambda_{\text{hs}}h^{2}s^{2}\,.$
(1)

The positivity condition on the quartic part of the potential above is the
same as the $\mathbb{Z}_{2}$ symmetric potential and is given by
$\lambda_{\text{h}}>0$, $\lambda_{\text{h}}>0$ and

$\lambda_{\text{hs}}>0~{}\vee~{}\left(\lambda_{\text{hs}}<0~{}\wedge~{}\lambda^{2}_{\text{hs}}\leq\lambda_{\text{h}}\lambda_{\text{s}}\right)\,.$
(2)

The vacuum solutions for the potential in Eq. (1) can only have the structures
$(0,0)$, $(0,w)$ or $(v,w)$; the VEV solution $(v,0)$ is not allowed. Despite
the $\mathbb{Z}_{2}$ symmetric model ( i.e. the potential in Eq. (1) with
$\kappa_{\text{hs}}=\kappa_{\text{s}}=0$), all three extremum solutions can be
local minimum at the same time even if we take into account the positivity
condition. For the $\mathbb{Z}_{2}$ symmetric model after the electroweak
symmetry breaking the possible minima at the EW scale are either $(v,w)$ or
$(v,0)$, but for the generic potential in Eq. (1) the only vacuum solution at
the EW scale is inevitably $(v,w)$ which is given by,

$v=\frac{\sqrt{\mu^{2}_{\text{h}}-\kappa_{\text{hs}}w-\lambda_{\text{hs}}w^{2}}}{\sqrt{\lambda_{\text{h}}}},\hskip
56.9055ptw=p+\xi_{-}^{1/3}+\xi+^{1/3},$ (3)

where

$\begin{split}&\xi_{\pm}=q\pm\sqrt{q^{2}+\left(r-p^{2}\right)^{3}},\\\
&q=\frac{b^{3}+9a(6\kappa_{\text{hs}}\mu^{2}_{\text{h}}a-bc)}{216a^{3}},~{}~{}~{}~{}~{}~{}p=\frac{b}{6a},~{}~{}~{}~{}~{}~{}~{}r=\frac{c}{6a}\,,\\\
&a=\lambda^{2}_{\text{hs}}-\lambda_{\text{h}}\lambda_{\text{s}},~{}~{}~{}~{}~{}b=2\kappa_{\text{s}}\lambda_{\text{h}}-3\kappa_{\text{hs}}\lambda_{\text{hs}},~{}~{}~{}~{}~{}c=\kappa^{2}_{\text{hs}}-2\lambda_{\text{hs}}\mu^{2}_{\text{h}}+2\lambda_{\text{h}}\mu^{2}_{\text{s}}\,.\end{split}$
(4)

From Eq. (4) real solutions for $v$ and $w$ requires $q^{2}>(p^{2}-r)^{1/3}$
for the reality of $w$ and
$\mu^{2}_{\text{h}}>0,\lambda_{\text{hs}}\leq-\kappa^{2}_{\text{hs}}/4\mu^{2}_{\text{h}}$
for the reality of $v$. The parameters $\mu^{2}_{\text{h}}$ and
$\mu^{2}_{\text{s}}$ can be fixed by the stationary conditions for $(v,w)$ at
a given scale $\mu=\Lambda$,

$\begin{split}&\mu^{2}_{\text{h}}=\lambda_{\text{h}}v^{2}+\lambda_{\text{hs}}w^{2}+\kappa_{\text{hs}}w\\\
&\mu^{2}_{\text{s}}=\lambda_{\text{s}}w^{2}+\lambda_{\text{hs}}v^{2}+\kappa_{\text{s}}w+\frac{\kappa_{\text{hs}}v^{2}}{2w}\end{split}$
(5)

where the Higgs VEV is fixed at $v=246$ GeV at the EW scale, and for the
scalar VEV we take benchmark values $w=100$ GeV and $w=1,10,100$ TeV at the EW
scale. Note that we choose the free parameters of the model to be
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}},\kappa_{\text{h}},\kappa_{\text{hs}}$;
we do not use a mixing angle as a new free parameter, instead we directly deal
with the various $(v,w)$ inputs from scratch and investigate the properties of
the model based on the chosen vacuum.

$w$ |  $100$ GeV |  $1$ TeV |  $10$ TeV |  $100$ TeV
---|---|---|---|---
$\lambda_{\text{h}}$ | $(0.02,0.24)$ | $(0.002,0.222)$ | $(0.00054,0.209)$ | $(8.5\times 10^{-5},0.096)$
$\lambda_{\text{s}}$ | $(0.002,0.97)$ | $(2.26\times 10^{-6},0.014)$ | $(2.54\times 10^{-6},0.00014)$ | $(3.10,11.45)\times 10^{-7}$
$\lambda_{\text{hs}}$ | $(-0.098,0.098)$ | $(-0.019,0.02)$ | $(-0.0021,0.0023)$ | $(-7.73,8.90)\times 10^{-6}$
$\kappa_{\text{s}}$ | $(-0.98,0.99)$ | $(-1,0.97)$ | $(-0.96,0.99)$ | $(-0.25,-0.05)$
$\kappa_{\text{hs}}$ | $(-1,1)$ | $(-0.95,0.99)$ | $(-0.96,0.98)$ | $(-0.71,0.95)$
$m_{s}$ |  $48$-$119$ GeV |  $46$-$118$ GeV |  $4$-$123$ GeV |  $8$-$109$ GeV
Table 2: The allowed region for the couplings
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}}$,
$\kappa_{\text{s}},\kappa_{\text{hs}}$ and $m_{s}$ for different singlet
scalar VEV benchmarks $w=0.1,1,10,100$ TeV, respecting the absolute vacuum
stability and positivity condition for the vacuum $(v,w)$ at the electroweak
scale with the assumptions $m_{H}=125$ GeV, $v_{H}=246$ GeV and $m_{s}<m_{H}$.

The Hessian matrix at a given scale $\mu=\Lambda$ in the scale-dependent
vacuum $(v,w)$ is given by,

$\mathcal{H}(v,w;\Lambda)=\left(\begin{matrix}3\lambda_{\text{h}}v^{2}+\lambda_{\text{hs}}w^{2}+\kappa_{\text{hs}}w-\mu^{2}_{\text{h}}&2\lambda_{\text{hs}}vw+\kappa_{\text{hs}}v\\\
2\lambda_{\text{hs}}vw+\kappa_{\text{hs}}v&3\lambda_{\text{s}}w^{2}+\lambda_{\text{hs}}v^{2}+2\kappa_{\text{s}}w-\mu^{2}_{\text{s}}\\\
\end{matrix}\right)$ (6)

where we have dropped in the matrix the scale-dependence of the couplings and
the VEVs. Taking into account the stationary condition on the vacuum $(v,w)$
in Eq. (5) the mass eigenvalues are,

$\begin{split}&m_{\pm}^{2}=v^{2}\lambda_{\text{h}}+w^{2}\lambda_{\text{s}}-\frac{v^{2}}{4w}\kappa_{\text{hs}}+\frac{w}{2}\kappa_{\text{s}}\pm\frac{v}{2}\times\\\
&\sqrt{\frac{v^{2}}{w^{2}}\left(\kappa_{\text{hs}}+4w\lambda_{\text{h}}\right)+\frac{w^{2}}{v^{2}}\left(\kappa_{\text{s}}+2w\lambda_{\text{s}}\right)+2\kappa_{\text{hs}}w\left(\lambda_{\text{s}}-8\lambda_{\text{hs}}\right)+4w\left(\kappa_{\text{s}}\lambda_{\text{h}}-4w\lambda_{\text{hs}}^{2}+2w\lambda_{\text{h}}\lambda_{\text{s}}\right)-4\kappa_{\text{hs}}^{2}}\\\
\end{split}$ (7)

where $m_{-}$ and $m_{+}$ can be either the Higgs or the scalar mass. In
section 4, we will consider both possibilities $m_{s}<m_{H}$ and $m_{s}>m_{H}$
at the EW scale. Although we set the initial inputs of the free parameters
such that $(v,w)$ is the absolute minimum at the EW scale, but running the
couplings in higher energy scales, the vacuum $(0,0)$ or $(0,w)$ may become
coexistent minima of the theory even deeper than the $(v,w)$, which results in
the instability of the vacuum at higher scales. To keep the $(v,w)$ to be the
absolute minimum at higher scales, we need to know the mass spectrum of other
possible vacua in higher scales. The mass matrix for the vacuum $(0,0)$ at a
given scale $\Lambda$ read

$\mathcal{M}(0,0;\Lambda)=\begin{pmatrix}-\mu^{2}_{\text{h}}&0\\\
0&-\mu^{2}_{\text{s}}\\\ \end{pmatrix}$ (8)

and the mass matrix for the vacuum $(0,w)$ is

$\mathcal{M}(0,w;\Lambda)=\begin{pmatrix}w^{2}\lambda_{\text{hs}}+w\kappa_{\text{hs}}-\mu^{2}_{\text{h}}&0\\\
0&2w^{2}\lambda_{\text{s}}+w\kappa_{\text{s}}\\\ \end{pmatrix}$ (9)

in which the stationary condition

$\mu^{2}_{\text{s}}=w^{2}\lambda_{\text{s}}+w\kappa_{\text{s}}$ (10)

for $(0,w)$ has been used. The initial values for the parameters
$\mu^{2}_{\text{h}}$ and $\mu^{2}_{\text{s}}$ in Eqs. (8) and (9) is given by
Eq. (5). In order for $(v,w)$ to stay the global minimum up to a desired scale
at least one of the mass eigenvalues in Eqs. (8) and (9) must be negative; in
this way $(v,w)$ would be an absolute minimum. Already at
$\Lambda=\mathcal{O}(m_{t})\sim 173$ GeV, from Eq. (8) we must require
$\mu^{2}_{\text{h}}>0$ or $\mu^{2}_{\text{s}}>0$ and from Eq. (9),
$\mu^{2}_{\text{h}}>w^{2}\lambda_{\text{hs}}+w\kappa_{\text{hs}}$ or
$2w^{2}\lambda_{\text{s}}+w\kappa_{\text{s}}<0$ in which both parameters
$\mu^{2}_{\text{s}}$ and $\mu^{2}_{\text{s}}$ are fixed at $\Lambda=173$ GeV
from Eq. (5).

##  3 Renormalizarion Group Equations

Figure 1: The viable region of input values for the couplings
$\lambda_{\text{h}}$, $\lambda_{\text{s}}$ and $\lambda_{\text{hs}}$ at the
electroweak scale for different singlet scalar VEV benchmarks $w=0.1,1,10,100$
TeV, respecting the absolute vacuum stability for the vacuum $(v,w)$, the
positinity condition and the perturbativity up to the Planck scale, with the
assumptions $m_{H}=125$ GeV, $v=246$ GeV and $m_{s}>m_{H}$.

The evolution of the couplings, fields and mass parameters in the model with
scale $\mu$ is given by the renormalization group equations (RGE). We extract
the RGEs for the model given in Eq. (1) up to one-loop using the Mathematica
package SARAH [21]. We take into account only the top quark Yukawa coupling
and ignore the couplings for the light quarks. The $\beta$-functions for the
gauge couplings and the Yukawa coupling are,

$\begin{split}16\pi^{2}&\beta_{g_{i}}=b_{i}g_{i}^{3}\\\
16\pi^{2}&\beta_{y_{t}}=\left(-\frac{17}{20}g^{2}_{1}-\frac{9}{4}g_{2}^{2}-8g_{3}^{2}\right)y_{t}+\frac{3}{2}y_{t}^{3}\\\
\end{split}$ (11)

where $b_{1}=41/10,~{}b_{2}=-19/6,~{}b_{3}=-7$. The $\beta$-functions
involving the Higgs and the singlet scalar couplings are given by,

$\begin{split}&16\pi^{2}\beta_{\lambda_{\text{h}}}=\left(\frac{27}{200}g_{1}^{4}+\frac{9}{20}g_{1}^{2}g_{2}^{2}+\frac{9}{8}g_{2}^{4}\right)+\left(-\frac{9}{5}g_{1}^{2}-9g_{2}^{2}+12y_{t}^{2}\right)\lambda_{\text{h}}+\frac{1}{2}\lambda_{\text{hs}}^{2}+24\lambda_{\text{h}}^{2}-6y_{t}^{4}\\\
&16\pi^{2}\beta_{\lambda_{\text{s}}}=2\lambda_{\text{hs}}^{2}+18\lambda^{2}_{\text{s}}\\\
&16\pi^{2}\beta_{\lambda_{\text{hs}}}=\left(-\frac{9}{10}g^{2}_{1}-\frac{9}{2}g_{2}^{2}+6\lambda_{\text{s}}+12\lambda_{\text{h}}+6y_{t}^{2}\right)\lambda_{\text{hs}}+4\lambda_{\text{hs}}^{2}\\\
&16\pi^{2}\beta_{\kappa_{\text{s}}}=6\kappa_{\text{hs}}\lambda_{\text{hs}}+18\lambda_{\text{h}}\kappa_{\text{s}}\\\
&16\pi^{2}\beta_{\kappa_{\text{hs}}}=\left(-\frac{9}{10}g_{1}^{2}-\frac{9}{2}g_{2}^{2}+4\lambda_{\text{hs}}+12\lambda_{\text{h}}+6y_{t}^{2}\right)\kappa_{\text{hs}}+2\kappa_{\text{s}}\lambda_{\text{hs}}\\\
\end{split}$ (12)

Figure 2: The plots show the running of the couplings
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}},\kappa_{\text{s}},\kappa_{\text{hs}}$
up to the Planck scale for different singlet scalar VEV benchmarks
$w=0.1,1,10,100$ TeV and with the constraints $m_{H}=125$ GeV, $v=246$ GeV and
$m_{s}>m_{h}$. The initial inputs for the couplings are chosen such that the
vacuum $(v,w)$ remains the absolute minimum respecting the positivity
condition and the perturbativity up to the Plack scale.

and the $\gamma$-functions for the VEVs and the mass parameters read,

$\begin{split}&16\pi^{2}\gamma_{v}=\left(\frac{9}{20}g_{1}^{2}+\frac{9}{4}g_{2}^{2}-3y_{t}^{2}\right)v\\\
&16\pi^{2}\gamma_{w}=0\\\
&16\pi^{2}\gamma_{\mu^{2}_{\text{h}}}=\left(-\frac{9}{10}g_{1}^{2}-\frac{9}{2}g_{2}^{2}+12\lambda_{\text{h}}+6y_{t}^{2}\right)\mu^{2}_{\text{h}}-2\kappa_{\text{hs}}^{2}+\lambda_{\text{hs}}\mu^{2}_{\text{s}}\\\
&16\pi^{2}\gamma_{\mu^{2}_{\text{s}}}=4\mu^{2}_{\text{h}}\lambda_{\text{hs}}+6\lambda_{\text{s}}\mu^{2}_{\text{s}}-4\kappa_{\text{hs}}^{2}-4\kappa_{\text{s}}^{2}\end{split}$
(13)

where the $\beta$-functions for a coupling $X$, and the $\gamma$-functions
(anomalous dimensions) for VEV or mass parameter $Y$, are defined as,

$\beta_{X}=\mu\frac{dX}{d\mu}\hskip
85.35826pt\gamma_{Y}=-\frac{\mu}{Y}\frac{dY}{d\mu}\,.$ (14)

Solving the RGEs in Eqs. (11), (12) and (13), the evolution of the couplings
and the VEVs will be known, so we can check the status of the stability, the
positivity and the perturbativity all together at any desired scale. In the
next section we will discuss the RGE solutions and the allowed parameters
inputs that we can use at the EW scale.

Figure 3: The viable region of input values for the couplings
$\lambda_{\text{h}}$, $\lambda_{\text{s}}$ and $\lambda_{\text{hs}}$ at the
electroweak scale for different singlet scalar VEV benchmarks $w=0.1,1,10,100$
TeV, respecting the absolute vacuum stability for the vacuum $(v,w)$, the
positinity condition and the perturbativity up to the Planck scale, with the
assumptions $m_{H}=125$ GeV, $v=246$ GeV and $m_{s}<m_{H}$.

##  4 Vacuum Stability, Positivity and Perturbativity

In this section we numerically solve the RG equations presented in section 3.
We always require that the vacuum $(v,w)$ defined in the EW scale
$\mathcal{O}(m_{t})\sim 173$ GeV, remains the absolute global minimum for
higher scales up to the Planck scale, hence pushing the stability up to the
Planck scale. Furthermore, we impose the positivity condition in Eq. (2)
(which is an scale-dependent condition evolving with the couplings) to hold
from the EW scale up to the Planck scale. We also discard the input values for
the parameters which lead to a Landau pole in the scales lower than the Planck
scale.

Figure 4: The plots show the running of the couplings
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}},\kappa_{\text{s}},\kappa_{\text{hs}}$
up to the Planck scale for different singlet scalar VEV benchmarks
$w=0.1,1,10,100$ TeV and with the constraints $m_{H}=125$ GeV, $v=246$ GeV and
$m_{s}<m_{h}$. The initial inputs for the couplings are chosen such that the
vacuum $(v,w)$ remains the absolute minimum respecting the positivity
condition and the perturbativity up to the Plack scale.

From the LHC experiments, the Higgs mass and the Higgs VEV are known: $v=246$
GeV and $m_{H}=125$ GeV. We will also fix the singlet scalar VEV by different
mass scale benchmarks $w=0.1,1,2,10,100$ TeV. In Eq. (7) one of the mass
eigenvalues is attributed to the Higgs mass. We investigate both cases
$m_{+}\equiv m_{s}>m_{H}\equiv m_{-}$ and $m_{-}\equiv m_{s}<m_{H}\equiv
m_{+}$. Among the parameters of the model as seen in section 2,
$\mu^{2}_{\text{h}}$ and $\mu^{2}_{\text{s}}$ are omitted by two stationary
conditions for $(v,w)$. Doing so we are left with the free independent
parameters being $\lambda_{\text{h}}$, $\lambda_{\text{s}}$,
$\lambda_{\text{hs}}$, $\kappa_{\text{s}}$ and $\kappa_{\text{hs}}$. The input
values for the free parameters at $\mathcal{O}(m_{t})$ scale should be chosen
such that the vacuum $(v,w)$ be the absolute global minimum. Any set of input
for the free parameters will give an input for $\mu^{2}_{\text{h}}$ and
$\mu^{2}_{\text{s}}$ from Eq. (5). Using these values in Eqs. (8) and (9) at
least one of the mass eigenvalues in $\mathcal{M}(0,0)$ and in
$\mathcal{M}(0,w)$ must be negative. Moreover, at the EW scale the input
values chosen for the free parameters should be bounded by the positivity
condition in Eq. (2), the positivity of the radicand in the mass expressions
in Eq. (13), and the positivity of the mass eigenvalue $m_{-}$ in Eq. (7).
Also depending on taking $m_{H}\equiv m_{+}\sim 125$ GeV or $m_{H}\equiv
m_{-}\sim 125$ GeV, the free parameters are bounded differently. The initial
values for the set of parameters satisfying the aforementioned constrained are
presented in Table 1 for $m_{s}>m_{H}$, and in Table 2 for $m_{s}<m_{H}$ with
the benchmarks $w=0.1,1,10,100$ TeV. In both tables the range of the allowed
singlet scale mass is shown. All the parameters $\lambda_{\text{h}}$,
$\lambda_{\text{s}}$, $\lambda_{\text{hs}}$, $\kappa_{\text{s}}$ and
$\kappa_{\text{hs}}$ are scanned in the interval $(-1,1)$.

After choosing a set of random input values for the parameters within the
regions in Tables 1 and 2, we solve numerically the RGEs given in section 3.
We repeat numeriously this process by taking input values and solving the RGEs
to cover all the allowed regions defined in Tables 1 and 2. Although the
initial values we found are suitable at scale $\Lambda\sim 173$ GeV, but as
running they vary in the higher scales and may violate one of the conditions
e.g. the stability constraint, the positivity condition or the perturbativity.
For a random set of inputs we check if all constraints are satisfied up to
Planck scale. In the case $m_{s}>m_{H}$ the viable initial values which lead
to an appropriate result is shown in Fig. 1 and the viable region of the
initial values in the case $m_{s}<m_{H}$ is shown in Fig. 3. The viable
regions given in Figs. 1 and (3) are for different benchmark singlet scalar
VEVs, $w=0.1,1,10,100$ TeV.

As seen in Fig. 1 for the case $w=100$ GeV, there is a narrow region which can
fulfill the desired condition up to Planck scale. But if we relax the Planck
scale implementm then there are a larger viable region for $w=100$ GeV. As we
increase the scalar VEV from $w=100$ GeV to $w=100$ TeV, we see in Fig. 1 that
the viable initial values for the couplings $\lambda_{\text{s}}$ and
$\lambda_{\text{s}}$ shrinks considerably. The coupling $\lambda_{\text{s}}$
remains $\mathcal{O}(1)$ for all singlet scalar VEV benchmarks.

In Fig. 3 that $m_{s}<m_{H}$ however the viable region for the case $w=100$
GeV is large and all coupling $\lambda_{\text{h}}$, $\lambda_{\text{s}}$ and
$\lambda_{\text{hs}}$ take $\mathcal{O}(1)$ initial values to fulfill the SPP
conditions up to Planck scale. As we increase the singlet scalar VEV from
$w=100$ GeV to $w=100$ TeV, the couplings $\lambda_{\text{s}}$ and
$\lambda_{\text{s}}$ must become smaller down to $\mathcal{O}(10^{-8})$ to
satisfy the SPP conditions. Also in Fig. 2 for $m_{s}>m_{H}$ and in Fig. 4 for
$m_{s}<m_{H}$ for each singlet scalar VEV benchmark $w$, we have also shown
the evolution of the free parameters for a randomly chosen set of parameters
within the regions in Figs 1 and 3. The singlet scalar mass is also bounded by
the SPP conditions. As seen in Tables 1 and 2 the singlet scalar mass varies
from about GeV up to about $0.5$ TeV.

##  5 Conclusion

The standard Model suffers from a vacuum metastability at high energy scale
around $10^{10}$ GeV. The addition of extra scalars in the hidden sector with
or without internal symmetries may stablize the vacuum up to the Planck scale.
However in general specially when symmetries are absent in the internal
configuration space (more investigations needed in future works), the
positivity condition might become strong enough to compete with the vacuum
stability at a given energy scale. As the simplest example we have
investigated a generic real singlet scalar extension of the Standard Model. We
have imposed the positivity condition for all scales alongside the absolute
vacuum stability and perturbativity to bound the free parameters of the model.
As seen in Tables 1 and 1, even before looking at higher scales, the free
parameters
$\lambda_{\text{h}},\lambda_{\text{s}},\lambda_{\text{hs}},\kappa_{\text{s}},\kappa_{\text{hs}}$
at the EW scale are strongly limited due respecting the absolute vacuum
stability and the positivity for the vacuum $(v,w)$, with $v=246$ GeV and
$w=0.1,1,10,100$ TeV being the Higgs and the singlet scalar vacuum expectation
values respectively. The bounds on the free parameter become stronger if we
want to keep the vacuum $(v,w)$ to be an absolute minimum, and at the same
time respecting the positivity condition and the perturbativity up to the
Planck scale (SPP conditions). The result for the viable space of input values
for the parameters to respect the SPP conditions is shown in Fig. 1 when
$m_{s}<m_{H}$ and in Fig. 3 when $m_{s}<m_{H}$. The upshot is that only for
$w=100$ and $m_{s}<m_{H}$ all the couplings are $\mathcal{O}(1)$. In other
cases, increasing the singlet scale VEV $w$ shrinks the viable region of the
couplings down to $\mathcal{O}(10^{-8})$ except for the coupling
$\lambda_{\text{h}}$ which remains in all case $\mathcal{O}(1)$ up to the
Planck scale.

We observe as well that the singlet scalar mass in the presence of the SPP
conditions can take values from $\mathcal{O}$(GeV) to $\mathcal{O}$(TeV). It
is worth to answer this question in future works that how the SPP conditions
act on the parameter space of a multi-scalar extension of the SM; whether the
SPP condition on multi-scalar theories are more or less restrictive than the
simple singlet scalar model we studied here.

## Acknowledgments

I would like to thank Alessandro Strumia for useful discussions. This work was
supported by the ERC grant NEO-NAT.

## References

  * [1] Georges Aad et al.  Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC.  Phys. Lett. B, 716:1–29, 2012.
  * [2] F. Abe et al.  Observation of top quark production in $\bar{p}p$ collisions.  Phys. Rev. Lett., 74:2626–2631, 1995.
  * [3] Fedor Bezrukov, Mikhail Yu. Kalmykov, Bernd A. Kniehl, and Mikhail Shaposhnikov.  Higgs Boson Mass and New Physics.  JHEP, 10:140, 2012.
  * [4] Dario Buttazzo, Giuseppe Degrassi, Pier Paolo Giardino, Gian F. Giudice, Filippo Sala, Alberto Salvio, and Alessandro Strumia.  Investigating the near-criticality of the Higgs boson.  JHEP, 12:089, 2013.
  * [5] N. Cabibbo, L. Maiani, G. Parisi, and R. Petronzio.  Bounds on the Fermions and Higgs Boson Masses in Grand Unified Theories.  Nucl. Phys. B, 158:295–305, 1979.
  * [6] Serguei Chatrchyan et al.  Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC.  Phys. Lett. B, 716:30–61, 2012.
  * [7] Chien-Yi Chen, S. Dawson, and I. M. Lewis.  Exploring resonant di-Higgs boson production in the Higgs singlet model.  Phys. Rev. D, 91(3):035015, 2015.
  * [8] Giuseppe Degrassi, Stefano Di Vita, Joan Elias-Miro, Jose R. Espinosa, Gian F. Giudice, Gino Isidori, and Alessandro Strumia.  Higgs mass and vacuum stability in the Standard Model at NNLO.  JHEP, 08:098, 2012.
  * [9] Joan Elias-Miro, Jose R. Espinosa, Gian F. Giudice, Hyun Min Lee, and Alessandro Strumia.  Stabilization of the Electroweak Vacuum by a Scalar Threshold Effect.  JHEP, 06:031, 2012.
  * [10] Yohei Ema, Mindaugas Karciauskas, Oleg Lebedev, Stanislav Rusak, and Marco Zatta.  Higgs–inflaton mixing and vacuum stability.  Phys. Lett. B, 789:373–377, 2019.
  * [11] Adam Falkowski, Christian Gross, and Oleg Lebedev.  A second Higgs from the Higgs portal.  JHEP, 05:057, 2015.
  * [12] P. M. Ferreira, R. Santos, and A. Barroso.  Stability of the tree-level vacuum in two Higgs doublet models against charge or CP spontaneous violation.  Phys. Lett. B, 603:219–229, 2004.  [Erratum: Phys.Lett.B 629, 114–114 (2005)].
  * [13] Emidio Gabrielli, Matti Heikinheimo, Kristjan Kannike, Antonio Racioppi, Martti Raidal, and Christian Spethmann.  Towards Completing the Standard Model: Vacuum Stability, EWSB and Dark Matter.  Phys. Rev. D, 89(1):015017, 2014.
  * [14] Karim Ghorbani and Parsa Hossein Ghorbani.  A Simultaneous Study of Dark Matter and Phase Transition: Two-Scalar Scenario.  JHEP, 12:077, 2019.
  * [15] Parsa Ghorbani.  Vacuum structure and electroweak phase transition in singlet scalar model.  10 2020.
  * [16] Parsa Ghorbani, Alessandro Strumia, and Daniele Teresi.  A landscape for the cosmological constant and the Higgs mass.  JHEP, 01:054, 2020.
  * [17] Matthew Gonderinger, Yingchuan Li, Hiren Patel, and Michael J. Ramsey-Musolf.  Vacuum Stability, Perturbativity, and Scalar Singlet Dark Matter.  JHEP, 01:053, 2010.
  * [18] Matthew Gonderinger, Hyungjun Lim, and Michael J. Ramsey-Musolf.  Complex Scalar Singlet Dark Matter: Vacuum Stability and Phenomenology.  Phys. Rev. D, 86:043511, 2012.
  * [19] Kristjan Kannike.  Vacuum Stability of a General Scalar Potential of a Few Fields.  Eur. Phys. J., C76(6):324, 2016.  [Erratum: Eur. Phys. J.C78,no.5,355(2018)].
  * [20] Valentin V. Khoze, Christopher McCabe, and Gunnar Ro.  Higgs vacuum stability from the dark matter portal.  JHEP, 08:026, 2014.
  * [21] Florian Staub.  Exploring new models in all detail with SARAH.  Adv. High Energy Phys., 2015:840780, 2015.