# Bayesian reconstruction of nuclear matter parameters from the equation of
state of neutron star matter

Sk Md Adil Imam1,2 N. K. Patra3 C. Mondal4 Tuhin Malik5 B. K. Agrawal1,2
<EMAIL_ADDRESS>1Saha Institute of Nuclear Physics, 1/AF Bidhannagar,
Kolkata 700064, India.  2Homi Bhabha National Institute, Anushakti Nagar,
Mumbai 400094, India.  3Department of Physics, BITS-Pilani, K. K. Birla Goa
Campus, Goa 403726, India  4Laboratoire de Physique Corpusculaire, CNRS,
ENSICAEN, UMR6534, Université de Caen Normandie, F-14000, Caen Cedex, France
5CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal

###### Abstract

The nuclear matter parameters (NMPs), those underlie in the construction of
the equation of state (EoS) of neutron star matter, are not directly
accessible. A Bayesian approach is applied to reconstruct the posterior
distributions of NMPs from EoS of neutron star matter. The constraints on
lower-order parameters as imposed by the finite nuclei observables are
incorporated through appropriately chosen prior distributions. The
calculations are performed with two sets of pseudo data on the EoS whose true
models are known. The median values of second or higher order NMPs show
sizeable deviations from their true values and associated uncertainties are
also larger. The sources of these uncertainties are identified as (i) the
correlations among various NMPs and (ii) leeway in the EoS of symmetric
nuclear matter, symmetry energy, and neutron-proton asymmetry which propagates
into the posterior distributions of the NMPs.

##  I Introduction

The bulk properties of neutron stars are instrumental in constraining the
equation of state (EoS) of dense matter Webb and Barret (2007); Lattimer and
Lim (2013). The conditions of charge neutrality and $\beta$-equilibrium
imposed on the neutron star matter renders it highly asymmetric leading to
neutron-proton ratio much larger than unity. The nuclear part of the EoS can
be decomposed into two main components: the EoS of symmetric nuclear matter
(SNM) and density-dependent symmetry energy. The knowledge of the EoS of
neutron star matter may provide an alternative probe to understand the
behavior of underlying symmetric nuclear matter and symmetry energy over a
wide range of density which may not be readily accessible in the terrestrial
laboratory. Usually, the components of neutron star matter EoS are expressed
in terms of nuclear matter parameters (NMPs), namely, the energy per nucleon
for symmetric nuclear matter, symmetry energy and their density derivatives
evaluated at the saturation density ($\rho_{0}\simeq 0.16$ fm-3). The lower-
order NMPs, governing the behavior of neutron star EoS at low densities are
determined by nuclear models calibrated to the bulk properties of finite
nuclei Lalazissis  _et al._ (1997); Piekarewicz (2006); Agrawal  _et al._
(2005); Klupfel  _et al._ (2009); Sulaksono  _et al._ (2009). The higher-order
NMPs are generally estimated using observed maximum neutron star mass together
with radius and tidal deformability corresponding to the neutron star with
canonical mass Zhang  _et al._ (2018); Cai and Li (2021); Gil  _et al._
(2021). Such investigations due to the lack of availability of enough
experimental data or for sake of simplicity are restricted to a small subspace
of NMPs.

Gravitational-wave astronomy through the observations of gravitational wave
signals emitted during the merging of binary neutron stars, promises
unprecedented constraints on the EoS of neutron star matter. The tidal
deformability inferred from these gravitational wave events encodes
information about the EoS. For the first time, gravitational wave event
GW170817 was observed by LIGO-Virgo detector from a low mass compact binary
neutron star (BNS) merger with the total mass of the system
$2.74_{-0.01}^{+0.04}M_{\odot}$ Abbott et al. (2019a, b). Another
gravitational wave signal likely originating from the coalescence of BNS
GW190425 is observed Abbott et al. (2020). These two events have already
triggered many theoretical investigations to constrain the EoS of neutron star
matter Abbott et al. (2017); Malik  _et al._ (2018); De  _et al._ (2018);
Fattoyev  _et al._ (2018); Landry and Essick (2019); Piekarewicz and Fattoyev
(2019); Malik  _et al._ (2019); Biswas  _et al._ (2020); Abbott et al. (2020);
Thi  _et al._ (2021a). The upcoming runs of LIGO-Virgo and the Einstein
Telescope are expected to observe many more gravitational wave signals emitted
from coalescing neutron stars. The mass and radius of neutron stars, observed
either in isolation or in binaries, by the Neutron star Interior Composition
Explorer Watts et al. (2016); Miller et al. (2019); Riley et al. (2019) have
offered complementary constraints on the EoS. A sufficiently large number of
such observations may be employed to constrain the NMPs directly which
underlie in the construction of the EoS of neutron star matter. Since one
needs to estimate simultaneously the values of about ten NMPs, investigations
along this direction require computationally efficient statistical tools which
allow the evaluation of the likelihood function for the experimental data that
may require appropriate marginalization.

A Bayesian approach is often applied to analyze gravitational-wave signals,
which involves nearly fifteen parameters, to infer their source properties
Ashton et al. (2019). It has been also extended to investigate the properties
of short gamma-ray burst Biscoveanu  _et al._ (2020), neutron star Coughlin
and Dietrich (2019); Hernandez Vivanco  _et al._ (2019); Biscoveanu  _et al._
(2019), the formation history of binary compact objects Lower  _et al._
(2018); Romero-Shaw  _et al._ (2019); Ramos-Buades  _et al._ (2020); Romero-
Shaw  _et al._ (2020); Zevin  _et al._ (2020), population using hierarchical
inference Abbott et al. (2019a); Talbot  _et al._ (2019) and to test general
relativity Keitel (2019); Ashton and Khan (2020); Payne  _et al._ (2019); Zhao
_et al._ (2019). Recently, Bayesian approach has become a useful statistical
tool for parameter estimation in the field of nuclear and nuclear-
astrophysics. It allows one to obtain joint posterior distributions of the
model parameters and the correlations among them for a given set of data.
Various constraints on the parameters known $\it{a}$ $\it{priori}$ are
incorporated through their prior distributions.

Extraction of nuclear matter properties from the chiral effective field theory
(EFT), in particular, the issue of over-fitting by appropriately choosing the
prior is described in great detail in Ref. Wesolowski  _et al._ (2016).
Recently, Bayesian techniques have also been employed to constrain symmetry
energy Somasundaram  _et al._ (2021), masses and radii of neutron stars
Drischler  _et al._ (2021) using the bounds obtained from chiral EFT. To
obtain the symmetry energy parameters from lower bound on neutron matter-
energy Tews  _et al._ (2017); to extract the crustal properties of neutron
star Carreau  _et al._ (2019); Thi  _et al._ (2021b); to limit the bounds on
cold neutron star matter EoS from observational constraints Riley  _et al._
(2018); Raaijmakers and Others (2020); Jiang  _et al._ (2020); to test the
compatibility of GW170817 event with multi-physics data Güven  _et al._
(2020); Biswas (2021a); to constrain neutron star matter from existing and
upcoming constraints on the gravitational wave and pulsar data Landry  _et
al._ (2020); to limit the neutron star EoS with microscopic and macroscopic
collisions Huth  _et al._ (2021); to filter models based on astrophysical
observations Biswas (2021b); or to limit the reach of nucleonic hypothesis in
the astrophysical context Thi  _et al._ (2021a); Mondal and Gulminelli (2021)
Bayesian techniques have been used extensively.

It is common to study different correlations in the posterior involving
parameters and observables alike, within a Bayesian analysis Biswas  _et al._
(2020); Güven  _et al._ (2020); Carreau (2020); Carreau  _et al._ (2019). The
origin of these uncertainties is embedded either in the underlying models used
or in the variances of the data employed. Existing correlations and
uncertainties among the extracted nuclear matter properties from a plethora of
nuclear physics data both from the laboratory as well as from the heavens and
theoretical calculations at low densities are often left for interpretation
Biswas  _et al._ (2020); Güven  _et al._ (2020); Carreau  _et al._ (2019). It
is also quite useful to test the limit of a certain type of data on physical
quantities which are extracted by employing Bayesian analyses. How far the
constraints on the static properties of a neutron star can pinpoint the
nuclear matter parameters is still a question of great interest. The same
applies to the data from heavy-ion collisions which can probe nuclear matter
at supra saturation densities Russotto  _et al._ (2011); Russotto and et. al.
(2016); Adamczewski-Musch and et. al. (2020). These studies can be done only
in a controlled environment, as the present observations associate large
uncertainties on the data. We have tried to do this by employing theoretical
modeling to mimic data on neutron star matter, as well as symmetric matter and
symmetry energy.

We first build the EoS of the neutron star matter by expanding it around
symmetric nuclear matter within the quadratic approximation as most commonly
employed. The EoS of symmetric nuclear matter and density-dependent symmetry
energy, the two main components, are further expanded around saturation
density within the Taylor and $\frac{n}{3}$ expansions. The expansion
coefficients in the former case are the individual NMPs and their linear
combinations in the latter case. A suitable set of NMPs is chosen so that the
resulting neutron star matter EoS is consistent with the currently observed
maximum mass of $\sim 2M_{\odot}$ and satisfies the causality condition. As
test cases, we employ these EoSs as pseudo data in a Bayesian analysis. The
true values of NMPs for the pseudo data are thus known $\it{a}$ $\it{priori}$.
The constraints imposed on the lower-order nuclear matter parameters by the
experimental data, for the bulk properties of finite nuclei, are incorporated
through appropriate choice of the prior distributions. The median values of
marginalized posterior distributions of NMPs and the associated uncertainties
on them as obtained for both the EoS models show similar trends. The inherent
nature of the model responsible for the deviations in the median values of
NMPs from their true values and uncertainties on them has been identified.

The paper is organized as follows, The Taylor and $\frac{n}{3}$ expansions for
the EoS of neutron star matter are briefly outlined in Sec. II. a Bayesian
approach is also discussed in the same section. The results for the posterior
distributions of NMPs obtained from the EoS of symmetric nuclear matter,
density-dependent symmetry energy and the EoS of neutron star matter is
presented in Sec. III. The main outcomes of the present investigation are
summarized in the last section.

##  II Methodology

The nuclear part of the energy per nucleon for neutron star matter
$\varepsilon(\rho,\delta)$ at a given total nucleon density, $\rho$ and
asymmetry, $\delta$ can be decomposed into the energy per nucleon for the SNM,
$\varepsilon(\rho,0)$ and the density-dependent symmetry energy, $J(\rho)$ in
the parabolic approximation as,

$\displaystyle\varepsilon(\rho,\delta)$ $\displaystyle=$
$\displaystyle\varepsilon(\rho,0)+J(\rho)\delta^{2}+...,$ (1)

where, $\delta=\left(\frac{\rho_{n}-\rho_{p}}{\rho}\right)$ with $\rho_{n}$
and $\rho_{p}$ being the neutron and proton densities, respectively. The value
of $\delta$ at a given $\rho$ is determined by the condition of
$\beta$-equilibrium and the charge neutrality. Once $\delta$ is known, the
fraction of neutron, proton, electron, muon can be easily evaluated. In the
following, we expand $\varepsilon(\rho,0)$ and $J(\rho)$ appearing in Eq. (1)
using Taylor and $\frac{n}{3}$ expansions. The coefficients of expansion in
case of the Taylor correspond to the individual nuclear matter parameters. In
the latter case, they are expressed as linear combinations of the nuclear
matter parameters. These EoSs are used as pseudo data in a Bayesian approach
to reconstruct the posterior distributions of nuclear matter parameters.

###  II.1 Taylor’s expansion

The $\varepsilon(\rho,0)$ and $J(\rho)$ can be expanded around the saturation
density $\rho_{0}$ as Chen  _et al._ (2005, 2009); Newton  _et al._ (2014);
Margueron  _et al._ (2018); Margueron and Gulminelli (2019),

$\displaystyle\varepsilon(\rho,0)$ $\displaystyle=$
$\displaystyle\sum_{n}\frac{a_{n}}{n!}\left(\frac{\rho-\rho_{0}}{3\rho_{0}}\right)^{n},$
(2) $\displaystyle J(\rho)$ $\displaystyle=$
$\displaystyle\sum_{n}\frac{b_{n}}{n!}\left(\frac{\rho-\rho_{0}}{3\rho_{0}}\right)^{n},$
(3)

so that,

$\displaystyle\varepsilon(\rho,\delta)$ $\displaystyle=$
$\displaystyle\sum_{n}\frac{1}{n!}(a_{n}+b_{n}\delta^{2})\left(\frac{\rho-\rho_{0}}{3\rho_{0}}\right)^{n},$
(4)

where the coefficients $a_{n}$ and $b_{n}$ are the nuclear matter parameters.
We truncate the sum in Eqs. (2) and (3) at 4th order, i.e., $n=0$ \- $4$.
Therefore, the coefficients $a_{n}$ and $b_{n}$ correspond to,

$\displaystyle a_{n}$ $\displaystyle\equiv$
$\displaystyle\varepsilon_{0},0,K_{0},Q_{0},Z_{0},$ (5) $\displaystyle b_{n}$
$\displaystyle\equiv$ $\displaystyle J_{0},L_{0},K_{\rm sym,0},Q_{\rm
sym,0},Z_{\rm sym,0}.$ (6)

In Eqs. (5) and (6), $\varepsilon_{0}$ is the binding energy per nucleon,
$K_{0}$ the incompressibility coefficient, $J_{0}$ the symmetry energy
coefficient, it’s slope parameter $L_{0}$, $K_{\rm sym,0}$ the symmetry energy
curvature parameter, $Q_{0}(Q_{\rm sym,0})$ and $Z_{0}(Z_{\rm sym,0})$ are
related to third and fourth order density derivatives of $\varepsilon(\rho,0)$
( $J(\rho)$), respectively. The subscript zero indicates that all the NMPs are
calculated at the saturation density.

It may be noticed from Eq. (4) that the coefficients $a_{n}$ and $b_{n}$ may
display some correlations among themselves provided the asymmetry parameter
depends weakly on the density. Further, the Eq. (4) may converge slowly at
high densities, i.e., $\rho>>4\rho_{0}$. This situation is encountered for the
heavier neutron stars. The neutron stars with a mass around $2M_{\odot}$,
typically have central densities $\sim 4-6\rho_{0}$.

###  II.2  ${\frac{n}{3}}$ expansion

Alternative expansion of $\varepsilon(\rho,\delta)$ can be obtained by
expanding $\varepsilon(\rho,0)$ and $J(\rho)$ as Lattimer and Prakash (2016);
Gil  _et al._ (2017),

$\displaystyle\varepsilon(\rho,0)$ $\displaystyle=$
$\displaystyle\sum_{n=2}^{6}(a^{\prime}_{n-2})\left(\frac{\rho}{\rho_{0}}\right)^{\frac{n}{3}},$
(7) $\displaystyle J(\rho)$ $\displaystyle=$
$\displaystyle\sum_{n=2}^{6}(b^{\prime}_{n-2})\left(\frac{\rho}{\rho_{0}}\right)^{\frac{n}{3}},$
(8) $\displaystyle\varepsilon(\rho,\delta)$ $\displaystyle=$
$\displaystyle\sum_{n=2}^{6}(a^{\prime}_{n-2}+b^{\prime}_{n-2}\delta^{2})\left(\frac{\rho}{\rho_{0}}\right)^{\frac{n}{3}}.$
(9)

We refer this as the ${\frac{n}{3}}$ expansion. It is now evident from Eqs.(7)
and (8) that the coefficients of expansion are no-longer the individual
nuclear matter parameters unlike in case of Taylor’s expansion. The values of
the NMPs can be expressed in terms of the expansion coefficients $a^{\prime}$
and $b^{\prime}$ as,

$\displaystyle\left(\begin{matrix}\varepsilon_{0}\\\ 0\\\ K_{0}\\\ Q_{0}\\\
Z_{0}\\\ \end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}1&1&1&1&1\\\ 2&3&4&5&6\\\ -2&0&4&10&18\\\
8&0&-8&-10&0\\\ -56&0&40&40&0\\\
\end{matrix}\right)\left(\begin{matrix}a^{\prime}_{0}\\\ a^{\prime}_{1}\\\
a^{\prime}_{2}\\\ a^{\prime}_{3}\\\ a^{\prime}_{4}\\\ \end{matrix}\right),$
(10) $\displaystyle\left(\begin{matrix}J_{0}\\\ L_{0}\\\ K_{\rm sym,0}\\\
Q_{\rm sym,0}\\\ Z_{\rm sym,0}\\\ \end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}1&1&1&1&1\\\ 2&3&4&5&6\\\ -2&0&4&10&18\\\
8&0&-8&-10&0\\\ -56&0&40&40&0\\\
\end{matrix}\right)\left(\begin{matrix}b^{\prime}_{0}\\\ b^{\prime}_{1}\\\
b^{\prime}_{2}\\\ b^{\prime}_{3}\\\ b^{\prime}_{4}\\\ \end{matrix}\right).$
(11)

The relations between the expansion coefficients and the NMPs are governed by
the nature of functional form for $\varepsilon(\rho,0)$ and $J(\rho)$. The
off-diagonal elements in the above matrices would vanish for the Taylor’s
expansion of $\varepsilon(\rho,0)$ and $J(\rho)$ as given by Eqs. (2) and (3),
respectively. Therefore, each of the expansion coefficients are simply the
individual NMPs given by Eqs. (5) and (6). Inverting the matrices in Eqs. (10)
and (11) we have,

$\displaystyle a^{\prime}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(360\varepsilon_{0}+20K_{0}+Z_{0}),$ $\displaystyle
a^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(-960\varepsilon_{0}-56K_{0}-4Q_{0}-4Z_{0}),$
$\displaystyle a^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(1080\varepsilon_{0}+60K_{0}+12Q_{0}+6Z_{0}),$
$\displaystyle a^{\prime}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(-576\varepsilon_{0}-32K_{0}-12Q_{0}-4Z_{0}),$
$\displaystyle a^{\prime}_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(120\varepsilon_{0}+8K_{0}+4Q_{0}+Z_{0}),$ (12)
$\displaystyle b^{\prime}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(360J_{0}-120L_{0}+20K_{\rm{sym,0}}+Z_{\rm{sym,0}}),$
$\displaystyle b^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{24}(-960J_{0}+328L_{0}-56K_{\rm{sym,0}}-4Q_{\rm{sym,0}}$
$\displaystyle-4Z_{\rm{sym,0}}),$ $\displaystyle b^{\prime}_{2}$
$\displaystyle=$
$\displaystyle\frac{1}{24}(1080J_{0}-360L_{0}+60K_{\rm{sym,0}}+12Q_{\rm{sym,0}}$
$\displaystyle+6Z_{\rm{sym,0}}),$ $\displaystyle b^{\prime}_{3}$
$\displaystyle=$
$\displaystyle\frac{1}{24}(-576J_{0}+192L_{0}-32K_{\rm{sym,0}}-12Q_{\rm{sym,0}}$
$\displaystyle-4Z_{\rm{sym,0}}),$ $\displaystyle b^{\prime}_{4}$
$\displaystyle=$
$\displaystyle\frac{1}{24}(120J_{0}-40L_{0}+8K_{\rm{sym,0}}+4Q_{\rm{sym,0}}$
(13) $\displaystyle+Z_{\rm{sym,0}}).$

Each of the coefficients $a^{\prime}$ and $b^{\prime}$ are the linear
combinations of nuclear matter parameters in such a way that the lower order
parameters may contribute dominantly at low densities. The effects of higher-
order parameters become prominent with the increase in density.

###  II.3 Bayesian estimation of NMPs

A Bayesian approach enables one to carry out detailed statistical analysis of
the parameters of a model for a given set of fit data. It yields joint
posterior distributions of model parameters which can be used to study not
only the distributions of given parameters but also to examine correlations
among model parameters. One can also incorporate prior knowledge of the model
parameters and various constraints on them through the prior distributions.
This approach is mainly based on the Bayes theorem which states that Gelman
_et al._ (2013),

$P(\bm{\theta}|D)=\frac{{\mathcal{L}}(D|\bm{\theta})P(\bm{\theta})}{\mathcal{Z}},$
(14)

where $\bm{\theta}$ and $D$ denote the set of model parameters and the fit
data. The $P(\bm{\theta}|D)$ is the joint posterior distribution of the
parameters, $\mathcal{L}(D|\bm{\theta})$ is the likelihood function,
$P(\bm{\theta})$ is the prior for the model parameters and $\mathcal{Z}$ is
the evidence. The posterior distribution of a given parameter can be obtained
by marginalizing $P(\bm{\theta}|D)$ over remaining parameters. The
marginalized posterior distribution for a parameter $\theta_{i}$ can be
obtained as,

$P(\theta_{i}|D)=\int P(\bm{\theta}|D)\prod_{k\not=i}d\theta_{k}.$ (15)

We use Gaussian likelihood function defined as,

$\displaystyle{\mathcal{L}}(D|\bm{\theta})$ $\displaystyle=$
$\displaystyle\prod_{j}\frac{1}{\sqrt{2\pi\sigma_{j}^{2}}}e^{-\frac{1}{2}\left(\frac{d_{j}-m_{j}(\bm{\theta)}}{\sigma_{j}}\right)^{2}}.$
(16)

Here the index $j$ runs over all the data, $d_{j}$ and $m_{j}$ are the data
and corresponding model values, respectively. The $\sigma_{j}$ are the adopted
uncertainties. The evidence $\mathcal{Z}$ in Eq. (14) is obtained by complete
marginalization of the likelihood function. It is relevant when employed to
compare different models. However in the present work $\mathcal{Z}$ is not
very relevant. To populate the posterior distribution of Eq. (14), we
implement a nested sampling algorithm by invoking the Pymultinest nested
sampling Buchner  _et al._ (2014) in the Bayesian Inference Library Ashton et
al. (2019).M

##  III Bayesian reconstruction of NMPs

We have considered Taylor and ${\frac{n}{3}}$ expansions in the previous
section to express the EoS for symmetric nuclear matter and the density-
dependent symmetry energy in terms of the NMPs. The EoS for neutron star
matter can thus be constructed for a given set of NMPs using Eqs. (4) and (9)
in a straightforward way. On the contrary, it is not evident that how reliably
the values of NMPs can be extracted once the EoS for the neutron star matter
is known. To illustrate, we construct EoS for the neutron star matter using
Taylor and ${\frac{n}{3}}$ expansion for a known set of NMPs. These EoSs are
then employed as pseudo data to reconstruct the marginalized posterior
distributions of the underlying NMPs through a Bayesian approach. Since the
true models for the pseudo data are known, the sources of uncertainties
associated with reconstructed NMPs may be analyzed more or less unambiguously.
A significant part of the uncertainties on the model parameters usually arises
from the intrinsic correlations among them Mondal  _et al._ (2015); Malik  _et
al._ (2018). The intrinsic correlations among the NMPs are the manifestation
of the various constraints imposed by the fit data Cai and Li (2021); Carreau
_et al._ (2019). These correlations may also depend on the choice of forms of
the functions for the EoS of symmetric nuclear matter and density-dependent
symmetry energy.

###  III.1  Likelihood function and prior distributions

To obtain the marginalized posterior distributions of model parameters within
a Bayesian approach one simply requires a set of fit data, a theoretical
model, and a set of priors for the model parameters as discussed in Sec. II.3.
The likelihood function for a given set of fit data is evaluated for a sample
of model parameters populated according to their prior distributions. The
joint posterior distributions of parameters are obtained with the help of the
product of the likelihood function and the prior distributions, Eq. (14). The
posterior distribution for individual parameters can be obtained by
marginalizing the joint posterior distribution with the remaining model
parameters. If the marginalized posterior distribution of a parameter is more
localized compared to its prior distribution, then, the parameter is said to
be well constrained by the fit data.

Table 1:  The values of nuclear matter parameters (in MeV) which are employed to construct various pseudo data using the Taylor and ${\frac{n}{3}}$ expansions. The parameters $\varepsilon_{0}$, $K_{0}$, $Q_{0}$ and $Z_{0}$ describes the EoS of the symmetric nuclear matter part and $J_{0}$, $L_{0}$, $K_{\rm sym,0}$, $Q_{\rm sym,0}$ and $Z_{\rm sym,0}$ describes density-dependent symmetry energy. The index ’N’ denotes the order of a given NMP.  N | Symmetric nuclear matter | Symmetry energy
---|---|---
0 | $\varepsilon_{0}$ | -16.0 | $J_{0}$ | 32.0
1 |  |  | $L_{0}$ | 50.0
2 | $K_{0}$ | 230 | $K_{\rm sym,0}$ | -100
3 | $Q_{0}$ | -400 | $Q_{\rm sym,0}$ | 550
4 | $Z_{0}$ | 1500 | $Z_{\rm sym,0}$ | -750

Our fit data are essentially the pseudo data for the EoS of symmetric nuclear
matter, density-dependent symmetry energy, and the EoS for neutron star matter
constructed from a suitable choice of NMPs as listed in Table 1. The values of
lower-order NMPs are close to those obtained from the SLy4 parameterization of
the Skyrme force calibrated to the bulk properties for a few selected finite
nuclei Chabanat  _et al._ (1997, 1998). The values of second or higher-order
NMPs are modified so that the EoS for the neutron star matter remains causal
for both the Taylor and $\frac{n}{3}$ expansions. Further, the maximum mass of
neutron stars for both the expansions satisfies the current lower bound of
$\sim 2M_{\odot}$. The likelihood function is obtained using Eq (16) for
pseudo data and the corresponding model values with the standard deviation,
$\sigma$ equal to unity at all densities ranging from 0.5 - 6$\rho_{0}$. The
present investigation may not be sensitive to the choice of the NMPs.

Table 2:  Two different sets P1 and P2 for the prior distributions of the nuclear matter parameters (in MeV). The parameters of Gaussian distribution (G) are $\mu$ (mean) and $\sigma$ (standard deviation). The parameters ’min’ and ’max’ denote the minimum and maximum values for the uniform distribution (U). The saturation density $\rho_{0}$ is taken to be 0.16 fm-3. Parameters | P1 | P2
---|---|---
| Pr-Dist | $\mu$ | $\sigma$ | Pr-Dist. | $\mu$ | $\sigma$
|  | min | max |  | min | max
$\varepsilon_{0}$ | G | -16 | 0.3 | G | -16 | 0.3
$K_{0}$ | G | 240 | 100 | G | 240 | 50
$Q_{0}$ | U | -2000 | 2000 | G | -400 | 400
$Z_{0}$ | U | -3000 | 3000 | U | -3000 | 3000
$J_{0}$ | G | 32 | 5 | G | 32 | 5
$L_{0}$ | U | 20 | 150 | G | 50 | 50
$K_{\rm sym,0}$ | U | -1000 | 1000 | G | -100 | 200
$Q_{\rm sym,0}$ | U | -2000 | 2000 | G | -550 | 400
$Z_{\rm sym,0}$ | U | -3000 | 3000 | U | -3000 | 3000

The calculations are performed for two different sets of priors. In Table 2,
we provide the details for the prior sets P1 and P2. Usually, if the
parameters are known only poorly, their prior distribution is taken to be
uniform. But, in case, if some information about a parameter is known $\it{a}$
$\it{priori}$, one simply assumes Gaussian distributions for the corresponding
parameter. The priors for $\varepsilon_{0}$, and $J_{0}$ are taken to be
Gaussian with their means and standard deviations consistent with the
constraints imposed by the finite nuclei properties. For most of the remaining
NMPs, the prior set P1 assumes uniform distributions. The prior set P2 further
imposes stronger constraints on the lower order parameters such as $K_{0}$ and
$L_{0}$ which are consistent with those obtained from finite nuclei
properties. The higher-order parameters are assumed to have wide Gaussian
distributions. We also impose an additional constraint on the symmetry energy
so that it always increases monotonically with density.

###  III.2 Symmetric nuclear matter and symmetry energy

The EoS of the symmetric nuclear matter $\varepsilon(\rho,0)$ and the density-
dependent symmetry energy $J(\rho)$ are the two main components which govern
the symmetric part and the deviations from it in the EoS of the neutron star
matter. The NMPs which are required in the constructions of
$\varepsilon(\rho,0)$ are $\varepsilon_{0},K_{0},Q_{0}$ and $Z_{0}$ and those
for $J(\rho)$ are $J_{0},L_{0},K_{\rm sym,0},Q_{\rm sym,0}$ and $Z_{\rm
sym,0}$. The nuclear matter parameters that appear in the expansions of
$\varepsilon(\rho,0)$ might be correlated with those NMPs appearing in the
expansion of $J(\rho)$ (cf Eqs. (4) and (9)). These correlations might prevent
the NMPs from being determined accurately. Moreover, the accurate values of
NMPs may also be masked by the strong correlations of symmetry energy with the
asymmetry parameter $\delta$ which determines the fractions of different
baryons and leptons at a given density. The sources of uncertainties in the
NMPs are intrinsically present in a EoS model. To avoid some of these
uncertainties, we first consider a Bayesian estimation of the NMPs for a given
$\varepsilon(\rho,0)$ and $J(\rho)$ separately, before embarking on their
estimations from a EoS of the neutron star matter.

Figure 1: Corner plots for marginalized posterior distributions of the NMPs
which appear in the expansion of the EoS for the symmetric nuclear matter
(top) and those in the density-dependent symmetry energy (bottom). The results
are obtained for the model M1 with the prior set P2. One dimensional posterior
distributions (light blue) plotted along the diagonal plots are also compared
with the corresponding prior distributions (light green). The vertical lines
indicate the 68% confidence interval of the NMPs. The confidence ellipses for
two-dimensional posterior distributions are plotted with $1\sigma$, $2\sigma$
and $3\sigma$ confidence intervals.

We perform a Bayesian analysis using Gaussian likelihood (cf. Eq. (16)) which
can be easily evaluated for a set of fit data together with corresponding
model values obtained for a sample of each parameter. We construct two sets of
pseudo data for $\varepsilon(\rho,0)$ and $J(\rho)$ each. These pseudo data
correspond to the Taylor and ${\frac{n}{3}}$ expansions referred hereafter as
models M1 and M2, respectively. The values of NMPs used for these pseudo data
are the same as listed in Table 1. So, the true values of NMPs for a given
pseudo data are known. The marginalized posterior distributions (PDs) for the
NMPs which underlie in the constructions of $\varepsilon(\rho,0)$ and
$J(\rho)$ are obtained separately.

The median values of NMPs and associated $1\sigma$ uncertainties from the
marginalized PDs as listed in Table 3 are obtained for the models M1 and M2
for two different prior sets. The median values of NMPs obtained for all the
different cases are very close to their true values as listed in Table 1. The
uncertainties on the NMPs obtained for both the prior sets are quite similar
to each other for a given model. However, the uncertainties are significantly
larger for the third and fourth-order NMPs in the case of model M2 in
comparison to the ones for M1. The uncertainties for third-order NMPs for M2
are about five times larger than those for the M1. It increases to more than
ten times for the fourth-order NMPs. The $1\sigma$ uncertainties obtained from
the marginalized PDs for NMPs as listed in Table 3 are significantly smaller
than the ones corresponding to their prior distributions. The prior
distribution for a given NMP thus appears relatively uniform compared to its
marginalized PD, as will be seen later.

Table 3: The median values and the $1\sigma$ errors for the nuclear matter parameters (in MeV) from their marginalized posterior distributions. The distributions of $\varepsilon_{0}$, $K_{0}$, $Q_{0}$ and $Z_{0}$ are reconstructed from the EoS of the symmetric nuclear matter and those for $J_{0}$, $L_{0}$, $K_{\rm sym,0}$, $Q_{\rm sym,0}$ and $Z_{\rm sym,0}$ from the density-dependent symmetry energy. The results are presented for the Taylor (M1) and ${\frac{n}{3}}$ (M2) expansions obtained using prior sets P1 and P2. NMPs | M1-P1 | M1-P2 | M2-P1 | M2-P2
---|---|---|---|---
$\varepsilon_{0}$ | $-16.0_{-0.2}^{+0.2}$ | $-16.0_{-0.2}^{+0.2}$ | $-16.0_{-0.2}^{+0.2}$ | $-16.0_{-0.2}^{+0.2}$
$K_{0}$ | $230_{-8}^{+8}$ | $230_{-8}^{+8}$ | $230_{-15}^{+14}$ | $230_{-14}^{+13}$
$Q_{0}$ | $-402_{-35}^{+35}$ | $-401_{-35}^{+35}$ | $-403_{-124}^{+128}$ | $-403_{-116}^{+122}$
$Z_{0}$ | $1502_{-53}^{+53}$ | $1501_{-53}^{+53}$ | $1515_{-773}^{+756}$ | $1517_{-739}^{+711}$
$J_{0}$ | $32.0_{-0.4}^{+0.4}$ | $32.0_{-0.4}^{+0.4}$ | $32.0_{-0.5}^{+0.5}$ | $32.0_{-0.4}^{+0.4}$
$L_{0}$ | $50.0_{-2.8}^{+2.6}$ | $50.0_{-2.6}^{+2.6}$ | $50.0_{-2.6}^{+2.6}$ | $50.0_{-2.5}^{+2.5}$
$K_{\rm sym,0}$ | $-100_{-18}^{+18}$ | $-100_{-18}^{+18}$ | $-100_{-27}^{+27}$ | $-100_{-24}^{+24}$
$Q_{\rm sym,0}$ | $551_{-59}^{+58}$ | $549_{-59}^{+58}$ | $548_{-193}^{+184}$ | $551_{-166}^{+163}$
$Z_{\rm sym,0}$ | $-750_{-75}^{+80}$ | $-749_{-77}^{+78}$ | $-734_{-1034}^{+1064}$ | $-759_{-906}^{+936}$

Figure 2: Same as Fig. 1, but, for model M2.

The marginalized PDs and the confidence ellipses for the NMPs, which determine
$\varepsilon(\rho,0)$ and $J(\rho)$, obtained for models M1 and M2 are
displayed as corner plots in Figs. 1 and 2. These results correspond to prior
set P2 which assumes very wide Gaussian distributions for the higher-order
NMPs (see Table 2). The one-dimensional marginalized PDs for the NMPs are
displayed along the diagonals of the corner plots (light blue full lines)
together with the corresponding prior distributions ( light green lines).
These PDs for all the NMPs are quite symmetric around their median values or
they represent Gaussian distribution. The PDs for most of the NMPs are
localized compared to the corresponding prior distributions. As a result, the
prior distributions for most of the NMPs appear to be flat in comparison to
the ones for the marginalized PDs. This is also the reflection of the
constraints imposed by the pseudo data. The median values of NMPs are very
close to their true values. The confidence ellipses are plotted along with the
off-diagonal elements of the corner plots corresponding to $1\sigma$,
$2\sigma$, and $3\sigma$ confidence intervals. The width and inclination of
the confidence ellipses for a pair of NMPs depend on their covariance which
determines the nature of the linear correlations among them Reinhard and
Nazarewicz (2010); Mondal  _et al._ (2015). It may be noted that the
correlation patterns obtained for both models are only marginally different.
However, the uncertainties in the higher-order parameters are significantly
larger for the model M2. This fact may be attributed to some complex intrinsic
correlations among the NMPs.It is clear from Eqs. (12) and (13) that the
expansion coefficients for the model M2 are the linear combinations of the
NMPs unlike those in M1. Furthermore, it can be seen that the higher-order
terms in model M2 relative to the lower-order ones have less impact as
compared to those in M1.

Figure 3:  The EoS for symmetric nuclear matter (top) and the symmetry energy
(bottom) as a function of density obtained within 95% confidence interval from
the posterior distributions of nuclear matter parameters for models M1 and M2
for the prior set P2. The pseudo data for M1 and M2 are shown by triangles and
squares, respectively.

In Fig.3, the variations of $\varepsilon(\rho,0)$ and $J(\rho)$ as a function
of density with $95\%$ confidence intervals are plotted. The $95\%$ confidence
interval lies in a very narrow range which once again points to the fact that
the large uncertainties on the NMPs are predominantly due to the correlations
among them. The results presented in Figs. 1-3 provides firm ground to perform
the analysis of the NMPs obtained in the following sub-section, using the EoS
of neutron star matter, which involves $\varepsilon(\rho,0)$ and $J(\rho)$,
simultaneously. On passing, we may also remark that though the values of
$\varepsilon(\rho,0)$ and $J(\rho)$ for the models, M1 and M2 are obtained
using the same set of NMPs, their behavior at high densities is significantly
different. The lower-order NMPs, which govern the low-density behavior of
$\varepsilon(\rho,0)$ and $J(\rho)$, maybe model-independent.

###  III.3 Neutron star matter

We now apply a Bayesian approach to reconstruct the marginalized PDs for the
NMPs using the EoS for the neutron star matter which satisfies the conditions
of $\beta$-equilibrium and charge neutrality. The EoS for neutron star matter
$\varepsilon(\rho,\delta)$ can be obtained using Eq. (1) for a given
$\varepsilon(\rho,0)$ and $J(\rho)$. We construct two sets of pseudo data for
$\varepsilon(\rho,\delta)$ corresponding to the models M1 and M2 obtained
using NMPs of Table 1 in Eqs. (4) and (9), respectively. The marginalized
distributions of all the nine NMPs are reconstructed simultaneously from the
pseudo data for $\varepsilon(\rho,\delta)$, since it consists of
$\varepsilon(\rho,0)$ and $J(\rho)$. A Bayesian analysis is performed with
models M1 and M2 for prior sets P1 and P2. The median values of the NMPs
obtained from the marginalized PDs and the corresponding $1\sigma$ errors are
listed in Table 4. The NMPs are somewhat better estimated for the prior set
P2. The symmetry energy slope parameter $L_{0}$ seems to be a special case as
its counterpart in the symmetric nuclear matter vanishes (cf. Eqs. (5) and
(6)). It may be noticed that the uncertainties on $K_{\rm sym,0}$, $Q_{\rm
sym,0}$ and $Z_{\rm sym,0}$ are much larger than their counterparts in the EoS
for symmetric nuclear matter. The errors on $K_{\rm sym,0}$ are, however,
similar to those derived from bulk properties of finite nuclei or other
correlation systematics Mondal  _et al._ (2018, 2017); Agrawal  _et al._
(2021), though, we have allowed larger variations of $K_{0}$ and $L_{0}$.
Once, the sufficiently accurate values of $\varepsilon(\rho,\delta)$ are
determined from various astrophysical observations, they can be combined with
finite nuclei constraints to obtain $L_{0}$ and $K_{\rm sym,0}$ in tighter
limits.

The results for NMPs in Table 4, obtained from the EoS of neutron star matter
are substantially different from those of Table 3, which were determined
separately from the EoS of symmetric nuclear matter and the density-dependent
symmetry energy. In general, these differences can be summarized as follows,
(i) the median values of NMPs in Table 4 show larger deviations from their
true values compared to those in Table 3, (ii) the uncertainties on the NMPs
determined from the EoS of neutron star matter are several times larger for
most of the NMPs, (iii) the uncertainties on the $Z_{0}$ and $Z_{\rm sym,0}$
in Table 4 are somewhat asymmetric about their median values reflecting their
non-Gaussian nature, and (iv) the ratios of uncertainties between the models,
M1 and M2 were obtained for third and fourth-order NMPs listed in Table 4 are
significantly smaller than the ones in Table 3. This already provides us some
clue that there are additional sources of uncertainties on the NMPs determined
from the EoS of the neutron star matter.

Table 4: Same as Table 3, but, the posterior distributions for all the nuclear matter parameters are reconstructed simultaneously from the EoS for the neutron star matter. NMPs | M1-P1 | M1-P2 | M2-P1 | M2-P2
---|---|---|---|---
$\varepsilon_{0}$ | $-16.0_{-0.3}^{+0.3}$ | $-16.0_{-0.3}^{+0.3}$ | $-16.0_{-0.3}^{+0.3}$ | $-16.0_{-0.3}^{+0.3}$
$K_{0}$ | $187_{-56}^{+65}$ | $221_{-28}^{+36}$ | $213_{-40}^{+47}$ | $230_{-25}^{+28}$
$Q_{0}$ | $-367_{-220}^{+196}$ | $-471_{-123}^{+113}$ | $-327_{-198}^{+243}$ | $-412_{-123}^{+159}$
$Z_{0}$ | $1518_{-236}^{+258}$ | $1632_{-157}^{+152}$ | $1307_{-1656}^{+1069}$ | $1637_{-1206}^{+835}$
$J_{0}$ | $31.8_{-2.6}^{+2.5}$ | $32.0_{-2.7}^{+2.6}$ | $32.0_{-2.5}^{+2.5}$ | $32.0_{-2.4}^{+2.6}$
$L_{0}$ | $52.8_{-19}^{+25}$ | $55.5_{-16}^{+17}$ | $53.1_{-19.3}^{+23.8}$ | $51.0_{-13.9}^{+14.0}$
$K_{\rm sym,0}$ | $-34_{-178}^{+142}$ | $-108_{-72}^{+76}$ | $-114_{-138}^{+113}$ | $-106_{-70}^{+70}$
$Q_{\rm sym,0}$ | $220_{-563}^{+755}$ | $486_{-264}^{+257}$ | $562_{-488}^{+572}$ | $522_{-241}^{+248}$
$Z_{\rm sym,0}$ | $807_{-1527}^{+1341}$ | $100_{-668}^{+876}$ | $-60_{-1921}^{+1944}$ | $-323_{-1643}^{+1920}$
Figure 4: Corner plots for the marginalized posterior distributions of nuclear
matter parameters (in MeV) obtained from the EoS for the neutron star matter
for the model M1 with prior set P2. The prior distributions (light green) are
also plotted for the comparison. Figure 5: Same as Fig. 4, but, for the model
M2.

It may be pointed out that there are several additional sources of
uncertainties on the NMIt is customary to study different correlations in the
posterior involving parameters and observables alike, within a Bayesian
analysis. …Ps which were avoided by reconstructing them separately from the
$\varepsilon(\rho,0)$ and $J(\rho)$ as shown in Figs. 1 and 2. These sources
of uncertainties are (i) inter-correlations of NMPs corresponding to
$\varepsilon(\rho,0)$ with the ones for $J(\rho)$, (ii) compensation in the
change of $J(\rho)$ with the asymmetry parameter $\delta$ and
$\varepsilon(\rho,0)$ in such a way that the EoS of neutron star matter
remains more or less unaltered. We analyze them in detail in the following.

The corner plots for the marginalized PDs for the NMPs in one and two
dimensions for the models M1 and M2 obtained for prior set P2 are displayed in
Figs. 4 and 5, respectively. The difference between the one-dimensional PDs
for the NMPs and corresponding prior distributions reflect the role of pseudo
data in constraining the NMPs. These marginalized posterior distributions of
the NMPs are at variance with those obtained separately from the EoS for the
symmetric nuclear matter and the density-dependent symmetry energy as shown in
Figs. 1 and 2. The shapes and the orientations of the confidence ellipses
suggest that the correlations among most of the pairs of NMPs have disappeared
or weakened. Strong correlations exist only between $K_{0}-Q_{0}$ ,
$Q_{0}-Z_{0}$ and $L_{0}-K_{\rm sym,0}$ pairs with correlation coefficient
$r\sim 0.8$ for model M1. However, in model M2 the $K_{0}-Q_{0}$ correlation
got disappeared. The inter-correlations of the NMPs corresponding to
$\varepsilon(\rho,0)$ with those for the $J(\rho)$ are almost absent. $Z_{0}$
and $Z_{\rm sym,0}$ show almost no correlation with the remaining NMPs.
Overall reduction in the correlations among the NMPs which are reconstructed
from the EoS of the neutron star matter, but, increase in their uncertainties
at the same time seem to be somewhat counter-intuitive. Other sources of
uncertainties as mentioned earlier need to be addressed.

Figure 6:  Plots of 68% and 95% confidence intervals for the EoS for neutron
star matter (top), the symmetric nuclear matter (middle), and the symmetry
energy (bottom) as a function of scaled density for model M2 with prior set
P2. The results are obtained from the posterior distributions of the NMPs
which are reconstructed from the pseudo data for the EoS of neutron star
matter( triangles). The spread in $\varepsilon(\rho,0)$ and $J(\rho)$ are
consistent with those for $\varepsilon(\rho,\delta)$.

We now examine the uncertainties in the NMPs which might arise due to the
allowed variations in the $\varepsilon(\rho,0)$, $J(\rho)$ and $\delta$ for a
given $\varepsilon(\rho,\delta)$. The value of asymmetry parameter $\delta$ is
mainly governed by the symmetry energy at a given density. As the symmetry
energy increases, the $\delta$ decreases. Thus, the symmetry energy and
$\delta$ may balance each other in such a way that the asymmetric part of the
EoS of neutron star matter remains unaffected. Moreover, the variations in the
asymmetric part of $\varepsilon(\rho,\delta)$ may also be compensated by the
symmetric nuclear matter $\varepsilon(\rho,0)$. In short, for a given
$\epsilon(\rho,\delta)$, the values of $J(\rho)$, $\varepsilon(\rho,0)$ and
$\delta$ may have some leeway. We use the marginalized PDs for the NMPs to
obtain 68% and 95% confidence intervals for $\varepsilon(\rho,\delta)$,
$\varepsilon(\rho,0)$ and $J(\rho)$ . The results are plotted only for the
M2-P2 case in Fig. 6. Other cases show similar qualitative trends and are not
shown here. The value of $\varepsilon(\rho,\delta)$ (top) vary in a narrow
bound at a given density, but, $\varepsilon(\rho,0)$ (middle) and $J(\rho)$
(bottom) have larger uncertainties. The $95\%$ confidence intervals for
$\varepsilon(\rho,0)$ and $J(\rho)$ are little asymmetric with respect to the
ones for the $68\%$ due to the non-Gaussian nature of higher-order NMPs as can
be seen from Table 4 and Fig. 5. The spread in the values of $J(\rho)$
increases with density rapidly beyond 2$\rho_{0}$. For $68\%$ confidence
interval, spread in $J(\rho)$ at 4$\rho_{0}$ is $\sim$ 36 MeV which increases
to $\sim$ 160 MeV at 6$\rho_{0}$, whereas, the spread in $\varepsilon(\rho,0)$
remains almost the same ($\sim$ 15 MeV) for the density in the range
4$\rho_{0}$ to 6$\rho_{0}$. The larger spread in $J(\rho)$ is balanced by
asymmetry parameter $\delta$ as well as by change in $\varepsilon(\rho,0)$
such that the EoS for neutron star matter remains almost unaffected.

Figure 7:  Plots of confidence ellipses with $1\sigma$ interval for the EoS
for symmetric nuclear matter (top) and asymmetry parameter (bottom) as a
function of symmetry energy at densities $\rho=4\rho_{0}$ and $6\rho_{0}$. The
symbol tilde denotes that the corresponding quantity is obtained with respect
to its median value.

These features are intrinsic in nature which is present in all the models for
neutron star matter (cf. Eq. 1). The marginalized PDs for the NMPs, plotted in
Figs. 4 and 5 effectively correspond to the values of $\varepsilon(\rho,0)$
and $J(\rho)$ displayed in Fig. 6. To probe further, the confidence ellipses
are plotted in Fig. 7 for the $68\%$ confidence intervals for
$\varepsilon(\rho,0)$, and $\delta$ as a function of $J(\rho)$ at fixed
densities $\rho=$ 4$\rho_{0}$ and 6$\rho_{0}$. The $J(\rho)$ is anti-
correlated with $\delta$ and $\varepsilon(\rho,0)$. The spread in the values
of $J(\rho)$ is predominantly due to its anti-correlation with $\delta$. The
uncertainties in $\varepsilon(\rho,0)$ and $J(\rho)$ propagate into the NMPs.
That is why, the marginalized posterior distributions of NMPs displayed in
Figs. 4 and 5 are significantly different in comparison to those shown in
Figs. 1 and 2. This also explains the reason behind the larger uncertainties
on the higher-order NMPs, which govern the high-density behavior of $J(\rho)$.
As the uncertainties in $\varepsilon(\rho,0)$ are smaller than those in
$J(\rho)$, they get reflected in the uncertainties of the corresponding NMPs
(see Table 4). It seems that the EoS of neutron star matter, usually
constrained by using several astrophysical observables alone may not be
sufficient to determine the NMPs in narrow bounds. The more reliable
determination would also require additional constraints on the EoS of
symmetric nuclear matter as well as on the density-dependent symmetry energy.
The experimental data on the EoS of symmetric nuclear matter from the heavy-
ion collision and the symmetry energy beyond the saturation density from the
isobaric analog states may help in constraining the NMPs further Danielewicz
_et al._ (2002); Tsang  _et al._ (2009); Danielewicz and Lee (2014); Horowitz
_et al._ (2014); Roca-Maza and Paar (2018).

Table 5:  Values of correlation coefficients among selected pairs of nuclear matter parameters obtained with the modified prior set P2′. The width of Gaussian priors for $Q_{0}$ and $L_{0}$ are reduced by factor of four and that for $K_{sym,0}$ by factor of two for P2′ in comparison to the prior set P2  Parameters | M1-P2 | M1-P2′ | M2-P2 | M2-P2′
---|---|---|---|---
$K_{0}-Q_{0}$ | -0.9 | -0.75 | -0.5 | -0.15
$Q_{0}-Z_{0}$ | \- 0.97 | -0.87 | -0.88 | -0.63
$K_{\rm sym,0}-L_{0}$ | -0.76 | -0.44 | -0.86 | -0.62
Table 6:  Same as 4, but, for a modified priors set P2′.  NMPs | M1-P2′ | M2-P2′ | NMPs | M1-P2′ | M2-P2′
---|---|---|---|---|---
$\varepsilon_{0}$ | $-16.0_{-0.3}^{+0.3}$ | $-16.0_{-0.3}^{+0.3}$ | $J_{0}$ | $32.0_{-2.0}^{+2.0}$ | $32.0_{-2.1}^{+2.2}$
|  |  | $L_{0}$ | $51.6_{-7.5}^{+7.4}$ | $50.0_{-6.9}^{+7.2}$
$K_{0}$ | $212_{-21}^{+24}$ | $225_{-21}^{+27}$ | $K_{\rm sym,0}$ | $-87_{-44}^{+45}$ | $-105_{-41}^{+43}$
$Q_{0}$ | $-425_{-55}^{+63}$ | $-408_{-59}^{+62}$ | $Q_{\rm sym,0}$ | $489_{-222}^{+219}$ | $544_{-209}^{+213}$
$Z_{0}$ | $1561_{-80}^{+75}$ | $1650_{-647}^{+573}$ | $Z_{\rm sym,0}$ | $-96_{-604}^{+775}$ | $-332_{-1621}^{+1857}$

We modify the prior distributions to simulate the influence of the constraints
on the NMPs derived from the data on the microscopic systems such as heavy-ion
collisions and the bulk properties of finite nuclei. These data are expected
to constraint the behavior of symmetric nuclear matter and symmetry energy
over a wide range of densities ranging from sub saturation density to supra
saturation densities up to 2-3$\rho_{0}$. The empirical values of pressure of
the symmetric nuclear matter at supra saturation densities may constrain the
value of $Q_{0}$. The data on iso-vector giant dipole resonance and neutron-
skin thickness in heavy nuclei may constraint the value of $L_{0}$ Roca-Maza
and Paar (2018). Once the values of $J_{0}$ and $L_{0}$ are constrained,
$K_{\rm sym,0}$ may also be somewhat constrained Tews  _et al._ (2017); Mondal
_et al._ (2017). We repeat our calculations by reducing the width of Gaussian
priors for $Q_{0}$, $L_{0}$ and $K_{\rm sym,0}$ in the prior set P2. For the
$Q_{0}$ and $L_{0}$, the values of width are reduced by a factor of four,
whereas, for the $K_{\rm sym,0}$ by a factor of two. In Table 5, we present
the values of correlation coefficients among some selected pairs of NMPs
obtained with the modified priors with those for the prior set P2. In general,
the correlations become weaker with the modified prior. Consequently, the
uncertainties on the NMPs have decreased as can be seen from Table 6. In
particular, the uncertainties on $Z_{0}$ has now become almost half for both
the models M1 and M2. The spread in the values of $\varepsilon(\rho,0)$ and
$J(\rho)$ become smaller by less than $10\%$ for the densities around
$2\rho_{0}$. But, their spreads at higher densities remain practically
unaltered. The issue presented in this paper needs further investigation. In
the present work, we have used the most commonly employed EoS expanded around
the symmetric nuclear matter. Some alternative representation of the EoS of
neutron star matter may be employed. One such form is the expansion of the EoS
around the neutron matter in powers of the proton fraction Forbes  _et al._
(2019).

##  IV Summary and outlook

A Bayesian approach has been applied to reconstruct the underlying nuclear
matter parameters which describe the EoS of the neutron star matter. The
calculations are performed using the EoS for neutron star matter by expanding
it around symmetric nuclear matter within the parabolic expansion as commonly
employed. The EoS of symmetric nuclear matter and density-dependent symmetry
energy required for such an EoS are expanded using Taylor and ${\frac{n}{3}}$
expansions. The expansion coefficients for the former are the individual
nuclear matter parameters and the linear combinations of them for the
${\frac{n}{3}}$ case. The pseudo data for the EoS for symmetric nuclear
matter, neutron star matter, and density-dependent symmetry energy are
constructed using both expansions. This pseudo data enable us to identify the
various sources of uncertainties associated with the marginalized posterior
distributions of NMPs, since, the true models are known. The posterior
distributions of the nuclear matter parameters are obtained using two
different sets of priors. One of the prior sets assumes that most of the
parameters are unknown, except for the lowest order ones which are the binding
energy per nucleon for the symmetric nuclear matter and the symmetry energy
coefficient at the saturation density.

The marginalized posterior distributions for the NMPs reconstructed separately
from the EoS of symmetric nuclear matter and density-dependent symmetry energy
are very much localized around their true values. But, the posterior
distributions for all the NMPs determined simultaneously from the EoS of
neutron star matter are at variance. The median values significantly deviate
from their true values and associated uncertainties are also larger, in
particular, for second or higher-order NMPs. The main sources of uncertainties
are found to be (i) the correlations among higher-order parameters describing
the EoS of symmetric nuclear matter and similar correlations in the case of
density-dependent symmetry energy and (ii) the larger uncertainties in the
symmetry energy at a given density due to its anti-correlation with asymmetry
parameter and the EoS of symmetric nuclear matter such that neutron star
matter EoS remains mostly unaffected. These are intrinsic in nature for the
EoS of neutron star matter obtained by expanding it around the symmetric
nuclear matter. The EoS of neutron star matter alone may not be sufficient to
determine the higher-order NMPs in narrow bounds. The higher-order NMPs are
correlated to the lower-order ones, thus, the low-density ab-initio
predictions for the EoS of symmetric nuclear matter and pure neutron matter
from the chiral effective field theory should also be considered for the
improved parameterizations. The experimental data on the EoS of symmetric
nuclear matter from the heavy-ion collision and the symmetry energy beyond the
saturation density from the isobaric analog states may further help in
constraining the nuclear matter parameters.

We have also performed the calculations by imposing stringent constraints on
prior distributions for $Q_{0}$, $L_{0}$ and $K_{sym,0}$ which led to the
reduction of the correlations among the NMPs. Consequently, the uncertainties
on some of the NMPs become smaller. The spread in the EoS of symmetric nuclear
matter and symmetry energy become smaller by less than $10\%$ for the
densities around $2\rho_{0}$. But, their spreads at higher densities remain
practically unaltered. It remains to be understood whether the sources of
various uncertainties identified in the present work are due to the expansion
of the EoS around the symmetric nuclear matter. It may be interesting to
perform an investigation using the EoS expanded around the pure neutron matter
instead of symmetric nuclear matterForbes  _et al._ (2019).

##  V Acknowledgements

The authors would like to thank Arunava Mukherjee, Saha Institute of Nuclear
Physics for useful discussion, a careful reading of the paper, and important
suggestions. N.K.P. would like to thank T.K. Jha for constant encouragement
and support and gratefully acknowledge the Department of Science and
Technology, India, for the support of DST/INSPIRE Fellowship/2019/IF190058.
C.M. acknowledges partial support from the IN2P3 Master Project “NewMAC”.

## References

  * Webb and Barret (2007) N. A. Webb and D. Barret, Astrophys. J. 671, 727 (2007).
  * Lattimer and Lim (2013) J. M. Lattimer and Y. Lim, Astrophys. J. 771, 51 (2013).
  * Lalazissis  _et al._ (1997) G. A. Lalazissis, J. König, and P. Ring, Phys. Rev. C 55, 540 (1997).
  * Piekarewicz (2006) J. Piekarewicz, Phys. Rev. Lett. 95, 122501 (2006).
  * Agrawal  _et al._ (2005) B. K. Agrawal, S. Shlomo, and V. K. Au, Phys. Rev. C 72, 014310 (2005).
  * Klupfel  _et al._ (2009) P. Klupfel, P. G. Reinhard, T. J. Burvenich, and J. A. Maruhn, Phys. Rev. C 79, 034310 (2009).
  * Sulaksono  _et al._ (2009) A. Sulaksono, T. J. Bürvenich, P. G. Reinhard, and J. A. Maruhn, Phys. Rev. C 79, 044306 (2009).
  * Zhang  _et al._ (2018) N.-B. Zhang, B.-A. Li, and J. Xu, Astrophys. J. 859, 90 (2018).
  * Cai and Li (2021) B. J. Cai and B. A. Li, Phys. Rev. C 103, 034607 (2021).
  * Gil  _et al._ (2021) H. Gil, P. Papakonstantinou, and C. H. Hyun, arXiv:2110.09802,[nucl-th] (2021).
  * Abbott et al. (2019a) B. P. Abbott et al., Phys. Rev. X 9, 011001 (2019a).
  * Abbott et al. (2019b) B. P. Abbott et al., Astrophys. J. 882, L24 (2019b).
  * Abbott et al. (2020) B. P. Abbott et al., Astrophys. J. 892, L3 (2020).
  * Abbott et al. (2017) B. P. Abbott et al., Phys. Rev. Lett. 119, 161101 (2017).
  * Malik  _et al._ (2018) T. Malik, N. Alam, M. Fortin, C. Providência, B. K. Agrawal, T. K. Jha, B. Kumar, and S. K. Patra, Phys. Rev. C 98, 035804 (2018).
  * De  _et al._ (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018).
  * Fattoyev  _et al._ (2018) F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019).
  * Piekarewicz and Fattoyev (2019) J. Piekarewicz and F. J. Fattoyev, Phys. Rev. C 99, 045802 (2019).
  * Malik  _et al._ (2019) T. Malik, B. K. Agrawal, J. N. De, S. K. Samaddar, C. Providência, C. Mondal, and T. K. Jha, Phys. Rev. C 99, 052801 (2019).
  * Biswas  _et al._ (2020) B. Biswas, P. Char, R. Nandi, and S. Bose, Phys. Rev. D 103, 103015 (2020).
  * Thi  _et al._ (2021a) H. D. Thi, C. Mondal, and F. Gulminelli, arXiv:2109.09675, [astro-ph.HE] (2021a).
  * Watts et al. (2016) A. L. Watts et al., Rev. Mod. Phys. 88, 021001 (2016).
  * Miller et al. (2019) M. C. Miller et al., Astrophys. J. 887, L24 (2019).
  * Riley et al. (2019) T. E. Riley et al., Astrophys. J. Lett. 887, L21 (2019).
  * Ashton et al. (2019) G. Ashton et al., Astrophys. J. Suppl. 241, 27 (2019).
  * Biscoveanu  _et al._ (2020) S. Biscoveanu, E. Thrane, and S. Vitale, Astrophys. J. 893, 38 (2020).
  * Coughlin and Dietrich (2019) M. W. Coughlin and T. Dietrich, Phys. Rev. D 100, 043011 (2019).
  * Hernandez Vivanco  _et al._ (2019) F. Hernandez Vivanco, R. Smith, E. Thrane, P. D. Lasky, C. Talbot, and V. Raymond, Phys. Rev. D 100, 103009 (2019).
  * Biscoveanu  _et al._ (2019) S. Biscoveanu, S. Vitale, and C.-J. Haster, Astrophys. J. 884, L32 (2019).
  * Lower  _et al._ (2018) M. E. Lower, E. Thrane, P. D. Lasky, and R. Smith, Phys. Rev. D 98, 083028 (2018).
  * Romero-Shaw  _et al._ (2019) I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Month. Not. R. Astron. Soc. 490, 5210 (2019).
  * Ramos-Buades  _et al._ (2020) A. Ramos-Buades, S. Husa, G. Pratten, H. Estellés, C. García-Quirós, M. Mateu-Lucena, M. Colleoni, and R. Jaume, Phys. Rev. D 101, 083015 (2020).
  * Romero-Shaw  _et al._ (2020) I. M. Romero-Shaw, N. Farrow, S. Stevenson, E. Thrane, and X. J. Zhu, Month. Not. R. Astron. Soc.: Lett. 496, L64 (2020).
  * Zevin  _et al._ (2020) M. Zevin, C. P. L. Berry, S. Coughlin, K. Chatziioannou, and S. Vitale, Astrophys. J. 899, L17 (2020).
  * Talbot  _et al._ (2019) C. Talbot, R. Smith, E. Thrane, and G. B. Poole, Phys. Rev. D 100, 043030 (2019).
  * Keitel (2019) D. Keitel, Res. Notes. AAS 3, 46 (2019).
  * Ashton and Khan (2020) G. Ashton and S. Khan, Phys. Rev. D 101, 064037 (2020).
  * Payne  _et al._ (2019) E. Payne, C. Talbot, and E. Thrane, Phys. Rev. D 100, 123017 (2019).
  * Zhao  _et al._ (2019) Z. C. Zhao, H. N. Lin, and Z. Chang, Chin. Phys. C 43, 075102 (2019).
  * Wesolowski  _et al._ (2016) S. Wesolowski, N. Klco, R. J. Furnstahl, D. R. Phillips, and A. Thapaliya, J. Phys. G 43, 074001 (2016).
  * Somasundaram  _et al._ (2021) R. Somasundaram, C. Drischler, I. Tews, and J. Margueron, Phys. Rev. C 103, 045803 (2021).
  * Drischler  _et al._ (2021) C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Phys. Rev. C 103, 045808 (2021).
  * Tews  _et al._ (2017) I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeitsev, Astrophys. J. 848, 105 (2017).
  * Carreau  _et al._ (2019) T. Carreau, F. Gulminelli, and J. Margueron, Phys. Rev. C 100, 055803 (2019).
  * Thi  _et al._ (2021b) H. D. Thi, T. Carreau, A. F. Fantina, and F. Gulminelli, Astron. Astrophys. 654, A114 (2021b).
  * Riley  _et al._ (2018) T. E. Riley, G. Raaijmakers, and A. L. Watts, Month. Not. R. Astron. Soc. 478, 1093 (2018).
  * Raaijmakers and Others (2020) G. Raaijmakers and Others, Astrophys. J. Lett. 893, L21 (2020).
  * Jiang  _et al._ (2020) J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, Astrophys. J. 892, 1 (2020).
  * Güven  _et al._ (2020) H. Güven, K. Bozkurt, E. Khan, and J. Margueron, Phys. Rev. C 102, 015805 (2020).
  * Biswas (2021a) B. Biswas, Astrophys. J. 921, 63 (2021a).
  * Landry  _et al._ (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020).
  * Huth  _et al._ (2021) S. Huth _et al._ , (2021), arXiv:2107.06229 [nucl-th] .
  * Biswas (2021b) B. Biswas, (2021b), arXiv:2106.02644 [astro-ph.HE] .
  * Mondal and Gulminelli (2021) C. Mondal and F. Gulminelli, (2021), arXiv:2111.04520 [nucl-th] .
  * Carreau (2020) T. Carreau,  _Modeling the (proto)neutron star crust : toward a controlled estimation of uncertainties_ , Ph.D. thesis (2020).
  * Russotto  _et al._ (2011) P. Russotto, P. Wu, M. Zoric, M. Chartier, Y. Leifels, R. Lemmon, Q. Li, J. Łukasik, A. Pagano, P. Pawłowski, and W. Trautmann, Physics Letters B 697, 471 (2011).
  * Russotto and et. al. (2016) P. Russotto and et. al., Phys. Rev. C 94, 034608 (2016).
  * Adamczewski-Musch and et. al. (2020) Adamczewski-Musch and et. al. (HADES Collaboration), Phys. Rev. Lett. 125, 262301 (2020).
  * Chen  _et al._ (2005) L.-W. Chen, C. M. Ko, and B.-A. Li, Phys. Rev. C 72, 064309 (2005).
  * Chen  _et al._ (2009) L.-W. Chen, B.-J. Cai, C. M. Ko, B.-A. Li, C. Shen, and J. Xu, Phys. Rev. C 80, 014322 (2009).
  * Newton  _et al._ (2014) W. G. Newton, J. Hooker, M. Gearheart, K. Murphy, D.-H. Wen, F. J. Fattoyev, and B.-A. Li, Eur. Phys. J. A 50, 41 (2014).
  * Margueron  _et al._ (2018) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018).
  * Margueron and Gulminelli (2019) J. Margueron and F. Gulminelli, Phys. Rev. C 99, 025806 (2019).
  * Lattimer and Prakash (2016) J. M. Lattimer and M. Prakash, Phys. Rept. 621, 127 (2016).
  * Gil  _et al._ (2017) H. Gil, P. Papakonstantinou, C. H. Hyun, T.-S. Park, and Y. Oh, Acta Phys. Polon. B 48, 305 (2017).
  * Gelman  _et al._ (2013) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, J. Carlin, H. Stern, D. Rubin, and D. Dunson,  _Bayesian Data Analysis Third edition_ , November (CRC Press, 2013).
  * Buchner  _et al._ (2014) J. Buchner, A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A. Merloni, M. Salvato, J. Donley, and D. Kocevski, Astron. Astrophys. 564, 015805 (2014).
  * Mondal  _et al._ (2015) C. Mondal, B. K. Agrawal, and J. N. De, Phys. Rev. C 92, 024302 (2015).
  * Chabanat  _et al._ (1997) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 627, 710 (1997).
  * Chabanat  _et al._ (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
  * Reinhard and Nazarewicz (2010) P.-G. Reinhard and W. Nazarewicz, Phys. Rev. C 81, 051303 (2010).
  * Mondal  _et al._ (2018) C. Mondal, B. K. Agrawal, J. N. De, and S. K. Samaddar, Int. J. Mod. Phys. E 27, 1850078 (2018).
  * Mondal  _et al._ (2017) C. Mondal, B. K. Agrawal, J. N. De, S. K. Samaddar, M. Centelles, and X. Viñas, Phys. Rev. C 96, 021302 (2017).
  * Agrawal  _et al._ (2021) B. K. Agrawal, T. Malik, J. N. De, and S. K. Samaddar, Eur. Phys. J. Spec. Top 230, 517 (2021).
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002).
  * Tsang  _et al._ (2009) M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and A. W. Steiner, Phys. Rev. Lett. 102, 122701 (2009).
  * Danielewicz and Lee (2014) P. Danielewicz and J. Lee, Nucl. Phys. A 922, 1 (2014).
  * Horowitz  _et al._ (2014) C. J. Horowitz, E. F. Brown, Y. Kim, W. G. Lynch, R. Michaels, A. Ono, J. Piekarewicz, M. B. Tsang, and H. H. Wolter, J. Phys. G 41, 093001 (2014).
  * Roca-Maza and Paar (2018) X. Roca-Maza and N. Paar, Prog. Part. Nucl. Phys. 101, 96 (2018).
  * Forbes  _et al._ (2019) M. M. N. Forbes, S. Bose, S. Reddy, D. Zhou, A. Mukherjee, and S. De, Phys. Rev. D 100, 083010 (2019).


# Single-junction quantum-circuit refrigerator

V. Vadimov  QCD Labs, QTF Centre of Excellence, Department of Applied Physics,
Aalto University, P.O. Box 13500, FI-00076 Aalto, Espoo, Finland  MSP Group,
QTF Centre of Excellence, Department of Applied Physics, Aalto University,
P.O. Box 11000, FI-00076 Aalto, Espoo, Finland  Institute for Physics of
Microstructures, Russian Academy of Sciences, 603950 Nizhny Novgorod, GSP-105,
Russia  A. Viitanen  QCD Labs, QTF Centre of Excellence, Department of Applied
Physics, Aalto University, P.O. Box 13500, FI-00076 Aalto, Espoo, Finland  T.
Mörstedt  QCD Labs, QTF Centre of Excellence, Department of Applied Physics,
Aalto University, P.O. Box 13500, FI-00076 Aalto, Espoo, Finland  T. Ala-
Nissila  MSP Group, QTF Centre of Excellence, Department of Applied Physics,
Aalto University, P.O. Box 11000, FI-00076 Aalto, Espoo, Finland
Interdisciplinary Centre for Mathematical Modelling, Department of
Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, UK  M.
Möttönen  QCD Labs, QTF Centre of Excellence, Department of Applied Physics,
Aalto University, P.O. Box 13500, FI-00076 Aalto, Espoo, Finland  VTT
Technical Research Centre of Finland Ltd., QTF Center of Excellence, P.O. Box
1000, FI-02044 VTT, Finland

###### Abstract

We propose a quantum-circuit refrigerator (QCR) based on photon-assisted
quasiparticle tunneling through a single normal-metal–insulator–superconductor
(NIS) junction. In contrast to previous works with multiple junctions and an
additional charge island for the QCR, we galvanically connect the NIS junction
to an inductively shunted electrode of a superconducting microwave resonator
making the device immune to low-frequency charge noise. At low characteristic
impedance of the resonator and parameters relevant to a recent experiment, we
observe that a semiclassical impedance model of the NIS junction reproduces
the bias voltage dependence of the QCR-induced damping rate and frequency
shift. For high characteristic impedances, we derive a Born–Markov master
equation and use it to observe significant non-linearities in the QCR-induced
dissipation and frequency shift. We further demonstrate that in this regime,
the QCR can be used to initialize the linear resonator into a non-thermal
state even in the absence of any microwave drive.

††preprint: AIP/123-QED

Superconducting quantum circuits [1, 2, 3, 4] have emerged as a highly
promising platform for quantum simulations and quantum information processing.
As these circuits advance in terms of performance and complexity, active
control of individual devices within the circuit constitutes a key challenge.
In particular, the control of temperature and dissipation in devices such as
transmission lines [5], microwave resonators [6] or qubits [7, 8, 9, 10, 11]
enables targeted cooling and initialization of these devices. This has already
been achieved with the invention of the quantum-circuit refrigerator (QCR)
[12], a stand-alone, on-chip device that can locally cool superconducting
circuits based on photon-assisted tunneling utilizing a pair of parallel
normal-metal–insulator–superconductor (NIS) junctions [13, 14]. We note that
the NIS junctions have been widely used for refrigeration of electronic
subsystem [15, 16, 17, 18, 19] due to the non-linear Peltier effect predicted
in Ref. 20. Contrary to this, the QCR cools down the photonic degrees of
freedom of the system.

In this work we present a simplified quantum-circuit refrigerator based on a
single NIS junction and analyze this device in the case of direct coupling to
a coplanar waveguide (CPW) resonator. The basic operation principle is similar
to that of a double-junction refrigerator: A small bias voltage
$V_{0}<\Delta/e$, where $\Delta$ is the superconductor gap parameter, is
applied across the junction. Quasiparticles can tunnel through the insulating
barrier, energetically allowed by absorption of photons from the coupled
circuit which decreases the temperature of the electric degrees of freedom of
the circuit [14, 21, 22, 23]. Such tunneling events induce quantum-state
transitions, for example, in a microwave resonator [24, 6]. With increasing
voltages $V_{0}>\Delta/e$, photon emission rate into the circuit approaches
that of photon absorption, which leads to heating [25, 26]. Therefore, we
operate the QCR in the range $\Delta-\hbar\omega_{\textrm{r}}<eV_{0}<\Delta$,
with $\omega_{\textrm{r}}$ being the angular frequency of the lowest mode of
the circuit.

We first present the design of a single-junction QCR coupled to a resonator,
compare it to conventional double-junction microcooler devices, and analyze
this system in terms of its tunable Lamb shift and dissipation rate. We
further investigate the case of a single-junction QCR coupled to a high-
impedance resonator, including multiphoton processes. Finally, we present our
conclusions and explore further theoretical and experimental opportunities
with this device.

The system under study in Fig. 1 consists of a single normal-
metal–insulator–superconductor tunnel junction galvanically coupled to a
superconducting resonator. The normal-metal lead is used to voltage bias the
junction, providing on-demand control of tunneling. This novel implementation
of a QCR differs from the previous devices [12] in two ways. First, we employ
a single NIS junction as opposed to a double junction in a SINIS structure.
Second, we avoid forming a charge island between the junction and the
resonator by galvanically coupling the resonator to the NIS junction as
opposed to capacitive coupling. The primary advantage of this approach is to
mitigate low-frequency charge noise caused by changes in the charge state of
the island and of the charge traps in its vicinity [27]. To achieve an equal
tunneling rate and cooling power as for the double-junction device [12, 14],
the corresponding single-junction QCR can be designed with half the tunneling
resistance. In an experimental setting, the single-junction QCR can be
controlled with a single line, whereas the double-junction QCR is typically
operated through two separate control lines, an input and an output line [12,
28, 6]. For the single-junction QCR, a galvanically coupled $\lambda/4$ CPW
resonator naturally fixes the superconducting electrode to the ground
potential at low frequencies.

Figure 1: (a) Design of a single-junction QCR galvanically coupled to a
$\lambda/4$ CPW resonator. The voltage profile of the fundamental mode at
angular frequency $\omega_{\mathrm{r}}$ has an anti-node at the superconductor
lead and a node where the center conductor connects to the ground plane
(green). The NIS junction of capacitance $C_{\mathrm{j}}$ is biased with a
voltage $V_{0}$ applied to the normal-metal lead. (b) Lumped-element circuit
model of the fundamental resonator mode coupled to the NIS junction.

We begin with an analysis of the dissipation and frequency shifts induced by
the coupling of the resonator to the NIS junction. To this end, we employ
classical equations of motion for the electromagnetic field in the circuit.
The current-voltage characteristic of an NIS junction is both strongly non-
linear and non-local in time. However, if we assume the ac component of the
voltage to be small, we can replace the junction by an effective frequency-
dependent conductance, which can be controlled by the dc component of the
voltage. To proceed, we define

$\displaystyle
V_{\textrm{NIS}}(t)=V_{\textrm{dc}}+\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}V_{\textrm{ac}}(\omega)\textrm{e}^{-\textrm{i}\omega
t}\;\mathrm{d}\omega,$ (1) $\displaystyle I_{\textrm{NIS}}(t)\approx
I_{\textrm{dc}}(V_{\textrm{dc}})+\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}G_{\textrm{NIS}}(V_{\textrm{dc}},\omega)V_{\textrm{ac}}(\omega)\textrm{e}^{-\textrm{i}\omega
t}\;\mathrm{d}\omega,$ (2)

where $G(V_{\textrm{dc}},\omega)$ is the conductance of the junction at bias
voltage $V_{\textrm{dc}}$ and angular frequency $\omega$. Such an expansion is
valid in the limit [29] $|eV_{\textrm{ac}}(\omega)|\ll|\hbar\omega|$.

In order to find the dissipation rate and the Lamb shift of the resonator
caused by the coupling to the NIS junction we write the Kirchhoff’s law for
the circuit shown in Fig. 1(b) as

$I_{\textrm{NIS}}[V_{0}-\dot{\varphi}]=C\ddot{\varphi}+\frac{\varphi}{L},$ (3)

where $C=C_{\textrm{j}}+C_{\textrm{r}}$,

$\varphi(t)=\int\limits_{-\infty}^{t}V(t^{\prime})\;\mathrm{d}t^{\prime},$ (4)

and $V(t)$ is the voltage across the capacitance $C_{\textrm{r}}$. Employing
Eq. (2) and going to the Fourier picture we obtain an equation for the
eigenfrequency of the resonator:

$-\omega^{2}C-\mathrm{i}\omega G_{\textrm{NIS}}(V_{0},\omega)+\frac{1}{L}=0$
(5)

Provided the conductance is small $G_{\textrm{NIS}}(V_{0},\omega)\ll
Z_{\textrm{r}}^{-1}$ and changes smoothly on the frequency scale of
$G_{\mathrm{NIS}}(V_{0},\omega_{\mathrm{r}})Z_{\textrm{r}}\omega_{\mathrm{r}}$,
where $Z_{\textrm{r}}=\sqrt{L/C}$ is the characteristic impedance of the
resonator, we can find the angular-frequency shift $\omega_{\textrm{L}}$ and
the dissipation rate $\gamma$ as

$\omega_{\textrm{L}}-\mathrm{i}\gamma\approx-\frac{\mathrm{i}}{2}{G_{\textrm{NIS}}(V_{0},\omega_{\textrm{r}})}\omega_{\textrm{r}}Z_{\textrm{r}},$
(6)

where $\omega_{\textrm{r}}=1/\sqrt{LC}$ is the resonance frequency of the bare
$LC$ circuit. Thus, both of these quantities are proportional to the imaginary
and real parts of the NIS junction conductance. We can use this simple
approach when the voltage across the capacitor satisfies the condition
$|eV(\omega)|\ll|\hbar\omega|$. In the quantum regime, voltage fluctuations
are bounded from below by zero point fluctuations [30] $\delta
V(\omega)/|\omega|\sim\sqrt{\hbar Z_{\textrm{r}}/2}$. Thus, we obtain the
applicability criterion for our result to be $\sqrt{{\pi
Z_{\textrm{r}}}/{R_{\textrm{K}}}}\ll 1$, where
$R_{\textrm{K}}=2\pi\hbar/e^{2}$ is the von Klitzing constant.

The real and the imaginary parts of the junction conductance are given by the
following expressions [29]:

$\mathop{\mathrm{Re}}\left[G_{\textrm{NIS}}(V_{0},\omega)\right]=\frac{1}{8R_{\textrm{j}}\omega}\int\limits_{-\infty}^{+\infty}\mathop{\mathrm{Im}}\left[\frac{\hbar\omega^{\prime}+\textrm{i}\gamma_{\textrm{D}}\Delta}{\sqrt{\Delta^{2}-(\hbar\omega^{\prime}+\textrm{i}\gamma_{\textrm{D}}\Delta)^{2}}}\right]\\\
\left\\{\tanh\left[\frac{\hbar(\omega^{\prime}+\omega)-eV_{0}}{2k_{\textrm{B}}T}\right]+\tanh\left[\frac{\hbar(\omega^{\prime}+\omega)+eV_{0}}{2k_{\textrm{B}}T}\right]-\right.\\\
\left.\tanh\left[\frac{\hbar(\omega^{\prime}-\omega)-eV_{0}}{2k_{\textrm{B}}T}\right]-\tanh\left[\frac{\hbar(\omega^{\prime}-\omega)+eV_{0}}{2k_{\textrm{B}}T}\right]\right\\}\;\mathrm{d}\omega^{\prime},$
(7)
$\mathop{\mathrm{Im}}\left[G_{\textrm{NIS}}(V_{0},\omega)\right]=\frac{1}{4R_{\textrm{j}}\omega}\int\limits_{-\infty}^{+\infty}\mathop{\mathrm{Re}}\left[\frac{\hbar\omega^{\prime}+\textrm{i}\gamma_{\textrm{D}}\Delta}{\sqrt{\Delta^{2}-(\hbar\omega^{\prime}+\textrm{i}\gamma_{\textrm{D}}\Delta)^{2}}}\right]\\\
\left(\frac{1}{2}\left\\{\tanh\left[\frac{\hbar(\omega^{\prime}+\omega)-eV_{0}}{2k_{\textrm{B}}T}\right]+\tanh\left[\frac{\hbar(\omega^{\prime}+\omega)+eV_{0}}{2k_{\textrm{B}}T}\right]+\right.\right.\\\
\left.\left.\tanh\left[\frac{\hbar(\omega^{\prime}-\omega)-eV_{0}}{2k_{\textrm{B}}T}\right]+\tanh\left[\frac{\hbar(\omega^{\prime}-\omega)+eV_{0}}{2k_{\textrm{B}}T}\right]\right\\}-\right.\\\
\left.\tanh\left[\frac{\hbar\omega^{\prime}-eV_{0}}{2k_{\textrm{B}}T}\right]-\tanh\left[\frac{\hbar\omega^{\prime}+eV_{0}}{2k_{\textrm{B}}T}\right]\right)\;\mathrm{d}\omega^{\prime},$
(8)

where $R_{\textrm{j}}$ is the junction resistance in the normal state,
$\gamma_{\textrm{D}}$ is the Dynes parameter of the junction, and $T$ is the
electron temperature of the normal metal.

Figure 2: (a) QCR-induced dissipation rate and (b) frequency shift of a CPW
resonator as functions of the bias voltage of the NIS junction. Blue lines
correspond to Eqs. (6), (7), and (8), the green dots show the experimental
data from Ref. 28. The QCR-induced frequency shift is defined to vanish at
zero voltage $V=0$. The device parameters are given by
$\omega_{\textrm{r}}=2\pi\times 4.67$ GHz, $Z_{\textrm{r}}=34.8$
$\mathrm{\Omega}$, $\Delta=0.215$ meV, $T=0.17$ K,
$\gamma_{\textrm{D}}=4\times 10^{-4}$, and the junction resistance
$R_{\textrm{j}}=16$ $\mathrm{k\Omega}$ is obtained from a fit to the
experimental relaxation rate. Panels (c) and (d) show the dissipation rate and
the frequency shift, respectively, for identical parameters except for $T=0$
K. The positions of the dips and peaks in the frequency shifts are determined
by $eV_{0}=\Delta,~{}\Delta\pm\hbar\omega_{\mathrm{r}}$.

The frequency shift and the dissipation rate of the resonator are shown in
Fig. 2 for parameters the corresponding to the experiments in Ref. [28] where
a double-junction QCR capacitively coupled to a $\lambda/2$ resonator is
employed. The theoretical results from Eqs. (6), (7), and (8) well agree with
the experimental data despite the differences in the design of the sample and
our model. Such a good correspondence can be explained by the fact that the
frequency shift and the dissipation rate are both proportional to the sum of
the conductances of the NIS junctions in the QCR. Provided that the SINIS
junction is symmetric, the expressions for the frequency shift and the
dissipation rate coincide with Eq. (6) where the tunneling resistance and the
characteristic impedance are adjusted accordingly.

If the condition $\sqrt{\pi Z_{\textrm{r}}/R_{\textrm{K}}}\ll 1$ is violated,
multi-photon processes become important in the quantum limit. In this case, we
to quantize the electromagnetic field in the circuit. Employing the standard
field quantization approach [24], we obtain the following Hamiltonian for the
circuit:

$\displaystyle\hat{H}=\hat{H}_{\textrm{r}}+\hat{H}_{\textrm{NS}}+\hat{H}_{\textrm{t}};$
(9)
$\displaystyle\hat{H}_{\textrm{r}}=\frac{\hat{q}^{2}}{2C}+\frac{\hat{\varphi}^{2}}{2L};$
(10)
$\hat{H}_{\textrm{NS}}=\sum\limits_{k\sigma}\left(\xi_{k}^{\textrm{n}}-eV_{0}\right)\hat{d}_{k\sigma}^{\dagger}\hat{d}_{k\sigma}+\sum\limits_{p\sigma}\xi_{p}^{\textrm{s}}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{p\sigma}+\\\
\Delta\sum\limits_{p}\left(\hat{c}_{p\uparrow}\hat{c}_{\bar{p}\downarrow}+\textrm{h.c.}\right);$
(11)
$\hat{H}_{\textrm{t}}=\sum\limits_{kp\sigma}\left(\gamma_{kp}\hat{d}_{k\sigma}^{\dagger}\hat{c}_{p\sigma}\textrm{e}^{\textrm{i}\frac{e\hat{\varphi}}{\hbar}}+\textrm{h.c.}\right),$
(12)

where $\hat{\varphi}$ and $\hat{q}$ are the canonically conjugate
$[\hat{\varphi},\hat{q}]=\textrm{i}\hbar$ flux and charge operators of the
$LC$ circuit, respectively, $\hat{d}_{k\sigma}$ is the annihilation operator
of an electron in mode $k$ with spin projection $\sigma$ in the normal metal,
$\xi_{k}^{\textrm{n}}$ is the energy of this mode with respect to the Fermi
energy, $\hat{c}_{p\sigma}$ is the annihilation operator of an electron in the
mode $p$ with spin projection $\sigma$ in the superconductor,
$\xi_{p}^{\textrm{s}}$ is the energy of this mode with respect to the Fermi
energy, and $\gamma_{kp}$ is the tunneling matrix element. We emphasize that
the coupling between the electromagnetic and fermionic degrees of freedom
occurs through the charge shift operator
$\exp(\textrm{i}{e\hat{\varphi}}/{\hbar})$ in the tunneling Hamiltonian, since
the tunneling of an electron through the junction is associated with a change
of the capacitor charge by a single elementary charge. Introducing the
notations

$\displaystyle\hat{a}=\frac{1}{\sqrt{2\hbar
Z_{\textrm{r}}}}\left(\hat{\varphi}+\textrm{i}Z_{\textrm{r}}\hat{q}\right);~{}\alpha=e\sqrt{\frac{Z_{\textrm{r}}}{2\hbar}}=\sqrt{\frac{\pi
Z_{\textrm{r}}}{R_{\textrm{K}}}};$ (13)
$\displaystyle\hat{F}=\textrm{e}^{\textrm{i}\alpha\left(\hat{a}+\hat{a}^{\dagger}\right)};~{}\hat{\Theta}=\sum\limits_{kp\sigma}\gamma_{kp}\hat{d}_{k\sigma}^{\dagger}\hat{c}_{p\sigma},$
(14)

we obtain a Redfield master equation [31] which governs the dynamics of the
resonator density operator $\hat{\rho}$ as

$\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}t}=\mathcal{L}(\hat{\rho})=-\frac{\textrm{i}}{\hbar}\left[\hat{H}_{\textrm{r}},\hat{\rho}\right]-\\\
\frac{1}{\hbar^{2}}\int\limits_{0}^{+\infty}\left\\{\left[\hat{F}\hat{F}^{\dagger}(-\tau)\hat{\rho}-\hat{F}^{\dagger}(-\tau)\hat{\rho}\hat{F}\right]\left\langle\hat{\Theta}^{\dagger}(\tau)\hat{\Theta}\right\rangle+\right.\\\
\left.\left[\hat{\rho}\hat{F}^{\dagger}(-\tau)\hat{F}-\hat{F}\hat{\rho}\hat{F}^{\dagger}(-\tau)\right]\left\langle\hat{\Theta}\hat{\Theta}^{\dagger}(\tau)\right\rangle+\right.\\\
\left.\left[\hat{F}^{\dagger}\hat{F}(-\tau)\hat{\rho}-\hat{F}(-\tau)\hat{\rho}\hat{F}^{\dagger}\right]\left\langle\hat{\Theta}(\tau)\hat{\Theta}^{\dagger}\right\rangle+\right.\\\
\left.\left[\rho\hat{F}(-\tau)\hat{F}^{\dagger}-\hat{F}^{\dagger}\hat{\rho}\hat{F}(-\tau)\right]\left\langle\hat{\Theta}^{\dagger}\hat{\Theta}(\tau)\right\rangle\right\\}\;\mathrm{d}\tau,$
(15)

where $\mathcal{L}$ is the Liouvillian of the system, and the time dependent
operators are given in the interaction picture:

$\displaystyle\hat{F}(\tau)=\textrm{e}^{\textrm{i}\frac{\hat{H}_{\textrm{r}}\tau}{\hbar}}\hat{F}\textrm{e}^{-\textrm{i}\frac{\hat{H}_{\textrm{r}}\tau}{\hbar}},~{}\hat{\Theta}(\tau)=\textrm{e}^{\textrm{i}\frac{\hat{H}_{\textrm{NS}}\tau}{\hbar}}\hat{\Theta}\,\textrm{e}^{-\textrm{i}\frac{\hat{H}_{\textrm{NS}}\tau}{\hbar}}.$
(16)

Quantum-mechanical averaging is taken over the biased Fermi–Dirac
distributions of the quasiparticles in the normal metal and in the
superconductor (see Supplementary Material (SM) for more information).

Instead of solving Eq. (15) we simplify it further by a so-called secular
approximation to the master equation, leaving only the terms which do not
oscillate in time in the interaction picture. Thus the Liouvillian
$\mathcal{L}$ can be replaced by its secular part $\mathcal{L}_{\textrm{sec}}$
(see SM for more information). The dynamics described by the original
Liouvillian $\mathcal{L}$ can violate the non-negativity of the density
operator and the secular approximation removes this inconvenience.
Furthermore, the secular Liouvillian $\mathcal{L}_{\textrm{sec}}$ has a block
structure, which significantly reduces the computational requirements. The
validity of both the Born–Markov and secular approximations is determined by
$\min(Z_{\textrm{r}},R_{\textrm{K}})/R_{\textrm{j}}\ll\min[1,k_{\textrm{B}}T/(\hbar\omega_{\textrm{r}})]$
(see SM for more information).

The steady-state density operator of the resonator satisfies
$\mathcal{L}_{\textrm{sec}}(\hat{\rho}_{0})=0$. Using the master equation we
calculate the response function $D(\omega)$ with respect to an infinitesimal
drive [32, 33, 34, 35] proportional to $\hat{a}+\hat{a}^{\dagger}$:

$D(\omega)=-\textrm{i}\mathop{\mathrm{Tr}}\left\\{\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\mathcal{L}_{\textrm{sec}}+\textrm{i}\omega\right)^{-1}\left[\hat{a}+\hat{a}^{\dagger},\hat{\rho}_{0}\right]\right\\}.$
(17)

For an isolated resonator in its ground state, the response function equals

$D_{0}(\omega)=\frac{2\omega_{\textrm{r}}}{\omega_{\textrm{r}}^{2}-(\omega+\textrm{i}0)^{2}}.$
(18)

The poles of the response function as functions of the complex-valued
frequency $\omega$ define the characteristic frequencies and decay rates of
the modes in the system. From the definition of $D(\omega)$, we observe that
these poles correspond to the eigenvalues $\lambda_{j}$ of the secular
Liouvillian $\mathcal{L}_{\textrm{sec}}$. Accordingly, the general response
function can be expanded as

$D(\omega)=\sum\limits_{j}\frac{f_{j}}{\lambda_{j}+\textrm{i}\omega},$ (19)

where the real and imaginary parts of $\lambda_{j}$ give the decay rate and
frequency of the $j^{\rm th}$ transition, respectively. The complex amplitude
$f_{j}$ quantifies how much this transition contributes to the response
function.

Figure 3: (a) QCR-induced dissipation rate, (b) frequency shift, and (c) the
mean photon number $\langle\hat{n}\rangle$ and the distance to the closest
thermal state $\eta$ for a high-impedance resonator as a function of the dc
voltage across the NIS junction. The parameters are given by
$\omega_{\textrm{r}}=2\pi\times 4.67$ GHz, $Z_{\textrm{r}}=20$
k$\mathrm{\Omega}$, $\Delta=0.215$ meV, $T=0.17$ K,
$\gamma_{\textrm{D}}=4\times 10^{-4}$, $R_{\textrm{j}}=640$
k$\mathrm{\Omega}$. (d)–(f) As panels (a)–(c) except for $T=0.02$ K. The color
of the lines in the panels (a), (b), (d), and (e) corresponds to the relative
absolute value of $f_{j}$ which quantifies the contribution of the resonance
to the response function $D(\omega)$.

Figure 3 shows the dissipation rate and the frequency shift of the resonator
as functions of the dc voltage applied to the NIS junction, obtained from
eigenvalues of the Liouvillian of the master equation. Here, we study a high-
impedance resonator with $Z_{\textrm{r}}=20$ k$\mathrm{\Omega}$, and
consequently increased the junction resistance to $R_{\textrm{j}}=640$
k$\mathrm{\Omega}$ in order to work in the validity regime of the Born–Markov
theory. In contrast to the low-impedance case, we observe here several
resonances that differ in frequency and decay rate. Lowering the temperature
to $T=0.02$ K significantly increases the number of pronounced branches. We
suggest that these branches correspond to the multi-photon transitions in the
resonator, similar to those found in Ref. 36. If the voltage is
$eV_{0}\approx\Delta-\hbar\omega_{\textrm{r}}n$, where $n$ is a positive
integer, then tunneling processes with absorption of $n$ and more photons from
the resonator are allowed. This implies that the resonator states with $n$ or
more photons rapidly decay to lower-energy states, while states with less than
$n$ photons decay relatively slowly. Thus, the system relaxes to some non-
equilibrium steady state with almost zero probability of having $n$ or more
photons, whereas the lower-energy states may have non-negligible occupation.
If $n=1$, the QCR absorbs all photons from the resonator and the system
relaxes toward the ground state. At negative $n$ (voltages above $\Delta/e$),
photon emission processes come into play and we effectively heat the system.
However, due to the complicated non-linear interaction between the resonator
and the QCR, the photon distribution remains non-thermal.

In order to confirm this qualitative picture, we calculate two static
properties of the steady state: the mean occupation number
$\langle\hat{n}\rangle=\mathop{\mathrm{Tr}}\\{\hat{\rho}_{0}\hat{a}^{\dagger}\hat{a}\\}$
and trace distance to the closest thermal state [37]
$\eta=\inf_{\beta}\left\|\hat{\rho}_{0}-\left(1-e^{-\beta}\right)e^{-\beta\hat{a}^{\dagger}\hat{a}}\right\|$,
where $\|\hat{A}\|=\mathop{\mathrm{Tr}}\\{\sqrt{\hat{A}\hat{A}^{\dagger}}\\}$.
These quantities are shown in Figs. 3(c) and 3(f) as functions of the bias
voltage $V_{0}$. At very low voltages the system is in a thermal state; hence
$\eta$ is very low and $\langle\hat{n}\rangle$ corresponds to the thermal
distribution with temperature $T$. With increasing voltage the resonator heats
up due to the finite Dynes parameter of the superconductor and weak tunneling
to the subgap states. When approaching the voltage $eV_{0}\sim\Delta$, the
occupation number decreases and so does the the distance to the thermal state,
since single-photon processes come into play. This effect is especially
pronounced in the low-temperature case shown in Fig. 3(f), where a sudden drop
almost to the ground state in the interval
$\Delta-\hbar\omega_{\textrm{r}}<eV_{0}<\Delta+\hbar\omega_{\textrm{r}}$ is
visible. With a further voltage increase, photon emission comes into play and
the QCR starts to heat the resonator. Interestingly, the steady state in this
regime still remains non-thermal for high-characteristic-impedance resonators.

To summarize, in this paper we have proposed a simple design of a QCR
consisting of a single NIS junction galvanically coupled to the cooled system.
We showed that the QCR-induced dissipation rate and frequency shift of the
low-characteristic-impedance resonator can be found from a semiclassical
treatment of the dynamics of the electromagnetic field in the circuit. In
order to study the high-characteristic-impedance resonators and to obtain a
distribution of photons in the steady state, we employed Born–Markov master
equation approach to the resonator dynamics. We calculated the response
function, the poles of which give the characteristic frequencies and decay
rates of the resonator. We showed that a single bright resonance of the
response function for the low-characteristic-impedance resonators splits into
several resonances for the high-characteristic-impedance resonators, which
correspond to the transitions between the shifted consecutive levels of the
resonator. The short decay time of the high-energy levels due to the multi-
photon transitions present in this system leads to a significant broadening of
some of the transitions. We also calculated the mean photon number and the
distance between the steady state to the closest thermal state. We showed that
for the high-characteristic-impedance resonators the photon distribution
function is essentially non-thermal due to the multi-photon transitions. With
this single-junction device we can achieve all the desired functionalities of
a double-junction QCR. It can be used as a platform for studying open quantum
systems which exhibit nontrivial phenomena in the non-linear limit.

## Supplementary Material

See Supplementary Material for the correlation functions of the NIS junction
and the details of the secular approximation to the Redfield master equation.

###### Acknowledgements.

We thank Alexander Mel’nikov, Sergei Sharov, Dmitry Golubev, Gianluigi
Catelani, and Matti Silveri for useful discussions. This work was financially
supported by the European Research Council under Grant No. 681311 (QUESS),
Grant No. 957440 (SCAR), by the Academy of Finland under Grant No. 318937 and
under its Centres of Excellence Program (projects 312300 and 312298), and by
the Finnish Cultural Foundation. We acknowledge the computational resources
provided by the Aalto Science-IT project.

## References

  * Nakamura, Pashkin, and Tsai [1999] Y. Nakamura, Y. A. Pashkin, and J. Tsai, “Coherent control of macroscopic quantum states in a single-cooper-pair box,” Nature 398, 786–788 (1999).
  * Wallraff  _et al._ [2004] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, “Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics,” Nature 431, 162–167 (2004).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, “A quantum engineer’s guide to superconducting qubits,” Applied Physics Reviews 6, 021318 (2019), https://doi.org/10.1063/1.5089550 .
  * Blais  _et al._ [2021a] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, “Circuit quantum electrodynamics,” Rev. Mod. Phys. 93, 025005 (2021a).
  * Partanen  _et al._ [2016] M. Partanen, K. Y. Tan, J. Govenius, R. E. Lake, M. K. Mäkelä, T. Tanttu, and M. Möttönen, “Quantum-limited heat conduction over macroscopic distances,” Nature Physics 12, 460–464 (2016).
  * Sevriuk  _et al._ [2019] V. A. Sevriuk, K. Y. Tan, E. Hyyppä, M. Silveri, M. Partanen, M. Jenei, S. Masuda, J. Goetz, V. Vesterinen, L. Grönberg, and M. Möttönen, “Fast control of dissipation in a superconducting resonator,” Applied Physics Letters 115, 082601 (2019).
  * Hsu  _et al._ [2020] H. Hsu, M. Silveri, A. Gunyhó, J. Goetz, G. Catelani, and M. Möttönen, “Tunable refrigerator for nonlinear quantum electric circuits,” Physical Review B 101, 235422 (2020), publisher: American Physical Society.
  * Yoshioka and Tsai [2021] T. Yoshioka and J. S. Tsai, “Fast unconditional initialization for superconducting qubit and resonator using quantum-circuit refrigerator,” Applied Physics Letters 119, 124003 (2021), https://doi.org/10.1063/5.0057894 .
  * Basilewitsch  _et al._ [2021] D. Basilewitsch, J. Fischer, D. M. Reich, D. Sugny, and C. P. Koch, “Fundamental bounds on qubit reset,” Phys. Rev. Research 3, 013110 (2021).
  * Magnard  _et al._ [2018] P. Magnard, P. Kurpiers, B. Royer, T. Walter, J.-C. Besse, S. Gasparinetti, M. Pechal, J. Heinsoo, S. Storz, A. Blais, and A. Wallraff, “Fast and Unconditional All-Microwave Reset of a Superconducting Qubit,” Physical Review Letters 121, 060502 (2018).
  * Geerlings  _et al._ [2013] K. Geerlings, Z. Leghtas, I. M. Pop, S. Shankar, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, “Demonstrating a Driven Reset Protocol for a Superconducting Qubit,” Physical Review Letters 110, 120501 (2013).
  * Tan  _et al._ [2017] K. Y. Tan, M. Partanen, R. E. Lake, J. Govenius, S. Masuda, and M. Möttönen, “Quantum-circuit refrigerator,” Nature communications 8, 1–8 (2017).
  * Mörstedt  _et al._ [2021] T. F. Mörstedt, A. Viitanen, V. Vadimov, V. Sevriuk, M. Partanen, E. Hyyppä, G. Catelani, M. Silveri, K. Y. Tan, and M. Möttönen, “Recent developments in quantum-circuit refrigeration,” (2021), arXiv:2111.11234 [quant-ph] .
  * Silveri  _et al._ [2017] M. Silveri, H. Grabert, S. Masuda, K. Y. Tan, and M. Möttönen, “Theory of quantum-circuit refrigeration by photon-assisted electron tunneling,” Physical Review B 96, 094524 (2017).
  * Giazotto  _et al._ [2006a] F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, “Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications,” Rev. Mod. Phys. 78, 217–274 (2006a).
  * Kuzmin  _et al._ [2019] L. Kuzmin, A. Pankratov, A. Gordeeva, V. Zbrozhek, V. Shamporov, L. Revin, A. Blagodatkin, S. Masi, and P. de Bernardis, “Photon-noise-limited cold-electron bolometer based on strong electron self-cooling for high-performance cosmology missions,” Communications Physics 2, 1–8 (2019).
  * Nguyen  _et al._ [2013] H. Q. Nguyen, T. Aref, V. J. Kauppila, M. Meschke, C. B. Winkelmann, H. Courtois, and J. P. Pekola, “Trapping hot quasi-particles in a high-power superconducting electronic cooler,” New Journal of Physics 15, 085013 (2013).
  * Pekola  _et al._ [2004] J. Pekola, R. Schoelkopf, J. Ullom, _et al._ , “Cryogenics on a chip,” Physics Today 57, 41–47 (2004).
  * Muhonen, Meschke, and Pekola [2012] J. T. Muhonen, M. Meschke, and J. P. Pekola, Reports on Prgress in Physics 75, 046501 (2012).
  * Bardas and Averin [1995] A. Bardas and D. Averin, “Peltier effect in normal-metal–superconductor microcontacts,” Phys. Rev. B 52, 12873–12877 (1995).
  * Giazotto  _et al._ [2006b] F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, “Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications,” Reviews of Modern Physics 78, 217–274 (2006b).
  * Courtois  _et al._ [2014] H. Courtois, F. W. J. Hekking, H. Q. Nguyen, and C. B. Winkelmann, “Electronic coolers based on superconducting tunnel junctions: Fundamentals and applications,” Journal of Low Temperature Physics 175, 799–812 (2014).
  * Lowell  _et al._ [2013] P. J. Lowell, G. C. O’Neil, J. M. Underwood, and J. N. Ullom, “Macroscale refrigeration by nanoscale electron transport,” Applied Physics Letters 102, 082601 (2013), https://doi.org/10.1063/1.4793515 .
  * Ingold and Nazarov [1992] G.-L. Ingold and Y. V. Nazarov, “Charge tunneling rates in ultrasmall junctions,” in  _Single charge tunneling_ (Springer, 1992) pp. 21–107.
  * Masuda  _et al._ [2018] S. Masuda, K. Y. Tan, M. Partanen, R. E. Lake, J. Govenius, M. Silveri, H. Grabert, and M. Möttönen, “Observation of microwave absorption and emission from incoherent electron tunneling through a normal-metal–insulator–superconductor junction,” Scientific Reports 8, 3966 (2018).
  * Hyyppä  _et al._ [2019] E. Hyyppä, M. Jenei, S. Masuda, V. Sevriuk, K. Y. Tan, M. Silveri, J. Goetz, M. Partanen, R. E. Lake, L. Grönberg, and M. Möttönen, “Calibration of cryogenic amplification chains using normal-metal–insulator–superconductor junctions,” Applied Physics Letters 114, 192603 (2019).
  * Hsu  _et al._ [2021] H. Hsu, M. Silveri, V. Sevriuk, M. Möttönen, and G. Catelani, “Charge dynamics in quantum-circuit refrigeration: Thermalization and microwave gain,” AVS Quantum Science 3, 042001 (2021).
  * Silveri  _et al._ [2019] M. Silveri, S. Masuda, V. Sevriuk, K. Y. Tan, M. Jenei, E. Hyyppä, F. Hassler, M. Partanen, J. Goetz, R. E. Lake, _et al._ , “Broadband lamb shift in an engineered quantum system,” Nature Physics 15, 533–537 (2019).
  * Tucker and Feldman [1985] J. R. Tucker and M. J. Feldman, “Quantum detection at millimeter wavelengths,” Reviews of Modern Physics 57, 1055 (1985).
  * Blais  _et al._ [2021b] A. Blais, A. L. Grimsmo, S. Girvin, and A. Wallraff, “Circuit quantum electrodynamics,” Reviews of Modern Physics 93, 025005 (2021b).
  * Redfield [1965] A. Redfield, “The theory of relaxation processes,” in  _Advances in Magnetic and Optical Resonance_ , Vol. 1 (Elsevier, 1965) pp. 1–32.
  * Albert  _et al._ [2016] V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang, “Geometry and response of lindbladians,” Phys. Rev. X 6, 041031 (2016).
  * Ban [2015] M. Ban, “Linear response theory for open quantum systems within the framework of the abl formalism,” Quantum Studies: Mathematics and Foundations 2, 51–62 (2015).
  * Campos Venuti and Zanardi [2016] L. Campos Venuti and P. Zanardi, “Dynamical response theory for driven-dissipative quantum systems,” Phys. Rev. A 93, 032101 (2016).
  * Villegas-Martínez, Soto-Eguibar, and Moya-Cessa [2016] B. Villegas-Martínez, F. Soto-Eguibar, and H. Moya-Cessa, “Application of perturbation theory to a master equation,” Advances in Mathematical Physics 2016 (2016).
  * Viitanen  _et al._ [2021] A. Viitanen, M. Silveri, M. Jenei, V. Sevriuk, K. Y. Tan, M. Partanen, J. Goetz, L. Grönberg, V. Vadimov, V. Lahtinen, and M. Möttönen, “Photon-number-dependent effective lamb shift,” Physical Review Research 3, 033126 (2021).
  * Hillery [1987] M. Hillery, “Nonclassical distance in quantum optics,” Phys. Rev. A 35, 725–732 (1987).