# Jets in Common Envelopes: a low mass main sequence star in a red giant

Diego López-Cámara1, Fabio De Colle2, Enrique Moreno Méndez3,4, Sagiv Shiber5,
Roberto Iaconi6
1Cátedras CONACyT – Universidad Nacional Autónoma de México, Instituto de
Astronomía, AP 70-264, CDMX 04510, México
2Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, A.
P. 70-543 04510 D. F. México
3Facultad de Ciencias, Universidad Nacional Autónoma de México, A. P. 70-543,
CDMX 04510, México
4CCyT, UACM, Casa Libertad, Ermita Iztapalapa 4163, Lomas de Zaragoza,
Iztapalapa, 09620, CDMX, México
5Department of Physics and Astronomy, Louisiana State University, Baton Rouge,
LA, 70803 USA
6Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-
ku, Kyoto 606-8502, Japan E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present small-scale three-dimensional hydrodynamical simulations of the
evolution of a 0.3 $M_{\odot}$ main sequence star which launches two
perpendicular jets within the envelope of a 0.88 $M_{\odot}$ red giant. Based
on previous large-scale simulations, we study the dynamics of the jets either
when the secondary star is grazing, when it has plunged-in, or when it is
well-within the envelope of the red giant (in each stage for $\sim$11 days).
The dynamics of the jets through the common envelope (CE) depend on the
conditions of the environment as well as on their powering. In the grazing
stage and the commencement of the plunge self-regulated jets need higher
efficiencies to break out of the envelope of the RG. Deep inside the CE, on
the timescales simulated, jets are choked independently of whether they are
self-regulated or constantly powered. Jets able to break out of the envelope
of the RG in large-scale simulations, are choked in our small-scale
simulations. The accreted angular momentum onto the secondary star is not
large enough to form a disk. The mass accretion onto the MS star is 1-10% of
the Bondi-Hoyle-Littleton rate ($\sim 10^{-3}-10^{-1}$ M⊙ yr-1). High
luminosity emission, from X-rays to UV and optical, is expected if the jets
break out of the CE. Our simulations illustrate the need for inclusion of more
realistic accretion and jet models in the dynamical evolution of the CEs.

###### keywords:

Binaries: general – Binaries (including multiple): close Accretion, accretion
disks – Stars: jets Methods: numerical – Hydrodynamics –

††pubyear: 2020††pagerange: Jets in Common Envelopes: a low mass main sequence
star in a red giant–LABEL:lastpage

##  1 Introduction

The common envelope (CE) phase is a key stage during the evolution of a binary
system (Ivanova et al., 2013). It is a short-lived phase in which the envelope
of the more massive stellar component evolves into a Red Giant (RG) star,
fills its Roche-lobe, and engulfs the less massive star (Paczynski, 1976; Iben
& Livio, 1993). The less massive star, which may be a low mass main sequence
(MS) star, inspirals around the core of the most massive, reaches a small
orbital separation and may produce a merger, or lead to the ejection of the CE
and the emergence of a close binary system (De Marco & Izzard, 2017). Thus,
the CE is key in producing high-energy and transitory phenomena (HEAP) such as
type Ia supernovae (Chevalier, 2012; Postnov & Yungelson, 2014), short GRBs
(Vigna-Gómez, et al., 2020), long GRBs (Brown et al., 2007; Moreno Méndez et
al., 2011), and gravitational waves (Abbott et al., 2016), among others.

The CE evolution is dictated by the release of orbital energy as the secondary
orbits within the CE (the so-called “$\alpha$ formalism”, van den Heuvel,
1976) and by other energy sources (Ivanova et al., 2013). Ionization energy
(Han et al., 1995; Reichardt et al., 2020), long period pulsations (Clayton et
al., 2017), dust formation (Glanz & Perets, 2018; Iaconi et al., 2020),
accretion onto the secondary star (Chamandy et al., 2018), and the launching
of jets (Soker, 2004) may have an important role in the CE evolution.

Angular momentum ($J$) is an important ingredient in stellar evolution. In
binary systems, the mass transferred from one star to another carries $J$ and
thus forms an accretion disk (around the accreting stellar component), from
which jets may be launched due to magnetic processes such as the Blandford &
Payne (1982) mechanism. Jets may shape planetary nebulae (see Soker, 2020, and
references therein), and bipolar nebula may be produced (Kamiński et al.,
2021). In isolated stars jets may play an important role in the explosions of
core collapse supernovae (Khokhlov et al., 1999; Papish & Soker, 2011, 2014;
Gilkis & Soker, 2014; Gilkis, Soker, & Papish, 2016; Bear, Grichener, & Soker,
2017; Piran et al., 2017), or may produce the ears observed in some core
collapse supernova remnants (Castelletti et al., 2006; González-Casanova et
al., 2014; Grichener & Soker, 2017), to mention a few examples.

If the jets are launched when the secondary is in the outskirts of the RG,
then the jets will graze the envelope of the more massive progenitor. The jets
may remove the envelope mass efficiently so that no CE would be produced and a
grazing envelope evolution would take place instead (GEE, see Soker, 2015;
Shiber, Kashi, & Soker, 2017). If the jets are launched within the envelope of
the RG, they may have an important role in the evolution of the CE phase
(Armitage & Livio, 2000; Soker, 2004; Chevalier, 2012; Shiber, Schreier &
Soker, 2016; Soker & Gilkis, 2018; Gilkis, Soker, & Kashi, 2019; Schreier,
Hillel & Soker, 2019; Soker, Grichener & Gilkis, 2019; Jones, 2020). The
Boomerang Nebula, for example, may have been produced by a CE event in which
jets played an important role (Sahai, Vlemmings, & Nyman, 2017).

Previous CE hydrodynamical (HD) simulations have mainly followed the large
scale evolution of the system, i.e., including the entire CE in the
computational domain (e.g., Passy et al., 2012; Kuruwita, Staff, & De Marco,
2016; Iaconi et al., 2017, 2019; Reichardt et al., 2019). Large scale
numerical simulation studies allow an understanding of the overall behavior of
the CE phase, the loss of material from the system and, eventually, the end of
the CE phase. Three-dimensional (3D) HD large scale studies have encountered
difficulties in unbinding large percentages of the stellar envelope and thus
to fully understand the end of the CE phase (Ricker & Taam, 2012; Ohlmann et
al., 2016; Staff et al., 2016; Chamandy et al., 2018, 2020; Prust & Chang,
2019; Sand et al., 2020; Glanz & Perets, 2021a, b). Also, large scale CE
studies involving the presence of jets show that a fraction of the envelope
may be blown away during the grazing envelope evolution (Shiber, Schreier &
Soker, 2016; Shiber & Soker, 2018), or when the MS star is immersed within the
envelope of the massive star (Shiber et al., 2019; Schreier, Hillel & Soker,
2019; Zou et al., 2022).

Small scale simulations, i.e. simulations with a computational box covering
only a fraction of the CE, are useful to understand the details of the
propagation of the jets, energy deposition, turbulence, accretion onto the
secondary, and the disk formation (e.g., MacLeod & Ramirez-Ruiz, 2015; MacLeod
et al., 2017; Chamandy et al., 2019; De et al., 2020; Everson et al., 2020).
3D HD small scale simulations of the effects of jets in the CE phase have
shown that accretion can be hypercritical and that the jets may play a key
role in the evolution of the system (Moreno Méndez, López-Cámara, & De Colle,
2017; López-Cámara, De Colle, & Moreno Méndez, 2019; López-Cámara, Moreno
Méndez, & De Colle, 2020). However, small scale simulations have the
disadvantage that only small time scales and the local behavior of the system
are followed ($~{}sim$ a few days, i.e. 10% of the time needed for the
companion star to plunge in the core of the primary; and $\sim 20\%$ of the
volume of the red giant in the simulations presented in this paper,
respectively).

In this paper we study, by a set of small scale high-resolution 3D HD
simulations, the propagation of jets launched from a MS star in the CE of a
RG. In this paper we zoom in one part of a simulation of Shiber et al. (2019,
S19 from now on) and study in detail the evolution of the jets and their
surroundings by taking advantage of the higher resolution as well as including
self-regulation in the jets (which has not been included in large scale CE
simulations). As initial conditions of the simulations, we use the CE large
scale configuration computed at different stages by S19. Three characteristic
stages are considered: a. when the MS is grazing the RG, b. when it is
starting to plunge-in through the envelope of the RG, and c. when the MS is
well within the CE. We consider constantly powered and self-regulated jets,
i.e., with a kinetic luminosity depending on the amount of matter accreting on
the secondary, as described below.

The paper is organised as follows. In Section 2 we describe the numerical
setups and the details of each model. In Section 3, we discuss the global
evolution, mass accretion and specific angular momentum rates during the
grazing and the plunge-in stages. In Section 4, we discuss our main results
and analyse the effects that successful or choked jets in the CE produce. In
Section 5, we present our conclusions.

##  2 Grazing and plunge-in setups and models

Figure 1 shows a schematic of the characteristic moments of the evolution of a
MS-RG system followed in this study. The MS star launches two perpendicular
jets as it spirals around the core of the RG. The phases which we study are
the grazing stage (panel a) and the plunge-in stage (panel b). For the latter,
we study the specific moments when the plunge-in has commenced, and when the
MS is well within the CE.The jets may be self-regulated by the accretion rate
that crosses the inner boundary ($R_{\rm in}$), or may be constantly powered.

Figure 1: Scheme of the characteristic moments of the evolution of a MS-RG
system followed in this study. Panel a) shows the grazing stage; and panel b)
the plunge-in stage, i.e. when the MS moves through the CE of the RG (we study
when the MS has just begun the plunge-in as well as when the MS is well within
the CE). We study self-regulated as well as constantly-powered jets, as well
as poorly- and highly-powered jets. They are launched perpendicular to the
equatorial plane from the inner boundary ($R_{\rm{in}}$, set at the orbital
separation $a_{\rm orb,init}$), and have an opening angle $\theta_{j}$ and
velocity $v_{j}$.

The initial configuration of the MS-RG system during the grazing, beginning of
the plunge-in (pib), and when the MS star is well within the plunge-in (piww)
stages were taken from S19. Also, the orbital separation between the MS and
the core of the RG ($a_{\rm orb}$) and hydro boundary conditions are set to
the S19 values at every time step. Specifically, we base our study on
simulation #5, in which the MS star orbits around the RG and does not launch
jets as it enters the CE 111Our simulations start with the companion at the
surface of the primary, however previous studies (Iaconi et al., 2017; Shiber
et al., 2019) showed that for a larger initial separation, the envelope is
more inflated at the point when the companion starts spiraling-in into the
envelope, and therefore we expect it would be easier for the jets to be
successful for larger initial separation. We chose simulation #5 of Shiber et
al. (2019) as a base to our simulations primarily because the orbital
evolution in simulations #5 (without jets) and #6 (with jets) is almost
identical in the first 50 days and thus it is easier to compare our results
with theirs.. The MS star has mass $m_{\rm MS}=0.30~{}\mbox{$M_{\odot}$}$, and
radius $R_{\rm MS}\sim 0.30~{}\mbox{$R_{\odot}$}$222The MS star is smaller
than the inner boundary ($R_{\rm MS}<R_{\rm{in}}$).. The RG has a total mass
of $m_{\rm RG}=0.88~{}\mbox{$M_{\odot}$}$, core mass $m_{\rm
RG,c}=0.39~{}\mbox{$M_{\odot}$}$, and initial radius of $R_{\rm
RG}=83~{}\mbox{$R_{\odot}$}=5.77\times 10^{12}$ cm. The initial orbital
separation of the grazing-envelope stage is $a_{\rm orb,init}\simeq R_{\rm
RG}$ (which corresponds to $t\equiv t_{0}=0$ days), at the pib stage is
$0.76R_{\rm RG}$ ($t=35$ days $\equiv t_{pib}$), and at the piww is
$0.48R_{\rm RG}$ ($t=52$ days $\equiv t_{piww}$). The evolution of the orbital
separation with time is shown in Figure 6 of S19. The change in the orbital
separation in each simulation is relatively small. We neglect the (small)
change in the orbital velocity caused by the movement of the frame of
reference at each time step.

Our simulations are performed in the rest frame of the MS star (non-inertial
frame), thus, a wind (due to the motion of the MS star around the RG), and the
Coriolis and Centrifugal forces are considered. The gravitational force is
estimated at each time step directly from the numerical simulation. The
gravitational effects of the RG star are taken as radial and consider the mass
enclosed within the orbital separation in the momentum and energy conservation
equations. The mass of the MS is set as as point mass in its correspondent
position (the amount of mass which it accretes during our integration time is
negligible compared to that of the MS star, thus the mass of the MS star is
maintained constant). We do not include self-gravity, but the interaction
between the gas and the companion star is followed precisely in the simulation
by S19, which includes self-gravity, and which our simulations are based upon.

In each of the MS-RG stages, we follow the evolution of self-regulated or
constantly powered jets, with different powering, during $10^{6}$ s (11.6
days). We assume that the material crossing the inner boundary is either
accreted onto the stellar surface or ejected as jets. The jets are injected
perpendicularly from an inner boundary located at $R_{\rm{in}}=10^{11}$ cm.
Each jet has an opening angle $\theta_{j}=30^{\circ}$, a velocity $v_{j}=438$
km s-1 and is defined inside a conical region with $R<R_{\rm{in}}$ and
$\theta<\theta_{j}$. The material inside the inner boundary is fully removed
at each time step (with density pressure and velocity rewritten to their
initial value within the inner boundary333As the gravitational force diverges
in a region which is outside the computational box, we do not include any
gravitational softening.).

The density of each jet is defined as
$\rho=\dot{M}/(4\pi(1-\cos{\theta_{j}})R_{\rm{in}}^{2}v_{j}$. The jet pressure
is defined as $P=nk_{B}T$, where $n=\rho/\mu m_{H}$ is the number density,
$\mu$ is the mean molecular weight, $m_{H}$ the proton mass, and we take a
temperature of $T/\mu=10^{6}$ K. The jets are injected without angular
momentum (although some angular momentum, as described in the paper, is
present at the inner boundary).

The jets were powered by a fraction of the mass accretion rate ($\eta_{\rm
sr}$ for the self-regulated jets and $\eta_{\rm c}$ for the constantly powered
jets). The luminosity of the self-regulating jets was $L_{j}=\eta_{\rm
sr}\dot{M}_{\rm acc}v_{j}^{2}$, where $\dot{M}_{\rm acc}$ is the mass
accretion rate that crosses $R_{\rm{in}}$. Meanwhile, the luminosity of the
constantly powered models was $L_{j}=\eta_{\rm c}\dot{M}_{\rm BHL}v_{j}^{2}$,
where $\dot{M}_{\rm BHL}$ is the assumed Bondi-Hoyle-Littleton rate (Hoyle &
Lyttleton, 1939; Bondi & Hoyle, 1944) at the RG surface in the simulations of
S19444The $\dot{M}_{\rm BHL}$ is estimated at the surface of the accretor, the
MS ($R\sim 0.3\mbox{$R_{\odot}$}$); the accreting wind velocity and density
are estimated from the Keplerian velocity (we assume there was no prior tidal
synchronization) and the density profile of the surface of the donor star,
when the MS star is near the surface of the RG star $\dot{M}_{\rm BHL}\approx
10^{25}$g s-1.. We explore accretion efficiency values up to 0.70 (the
high-$\eta$ values are used to explore the most extreme cases).

The simulations were carried out by employing the 3D, adaptive mesh
refinement, HD code Mezcal (De Colle et al., 2012). The numerical domain in
all stages covered a large fraction of the CE ($\Delta X=\Delta Y=0.76R_{\rm
RG}=4.42\times 10^{12}$ cm $,\Delta Z=2\Delta X$). Five levels of refinement,
corresponding to a maximum resolution of $\delta X=\delta Y=\delta
Z=4.90\times 10^{-2}R_{\rm RG}=3.44\times 10^{9}$ cm are used. A fixed
spherical region centered in the MS star and with a radius of $0.15\Delta Z$
is allowed to refine to the maximum resolution. The maximum levels of
refinement drop for regions outside $0.3\Delta Z,0.45\Delta Z,0.6\Delta Z$.
Each jet was injected in $\approx$ 700 cells at the solid angle produced by
$R=R_{\rm{in}}$ and $\theta\leq\theta_{j}$ in its respective hemisphere (the
reminder of the $R=R_{\rm{in}}$ surface, $\approx 10^{4}$cells, has an outlfow
condition). The inner part of the RG and the motion of the MS star were taken
into account by setting the $X_{\rm min}$ and the $Y_{\rm max}$ with the
values coming from the large scale simulation of S19. All other boundaries
were set with outflow boundary conditions.

The mass accretion rate ($\dot{M}$) and angular momentum rate ($\dot{J}$),
which cross $R=R_{\rm{in}}$, are saved before the material is removed from the
computational domain). The specific angular ($J_{\rm sp}$) momentum is
obtained by $J_{\rm sp}=\dot{J}$/$\dot{M}$. The integration time for each
simulation is t${}_{\rm int}=10^{6}$ s, and the Courant number is set at
$n_{\rm C}=0.1$. We verify that, in the absence of jets, the common envelope
remains in quasi-hydrostatic equilibrium for the duration of the simulation.
In the self-regulating jet models we have assumed that the material that
crosses $R_{\rm{in}}$ immediately alters the powering of the jets. The latter
is not true since the material requires a time-lag ($t_{\rm{lag}}$) in order
to reach the surface of the MS star. This however does not modify the main
results since including a $t_{\rm{lag}}>0$ does not affect the global outcome
(see for example, López-Cámara, De Colle, & Moreno Méndez, 2019).

Thirty models are followed, see Table 1. The models are labeled according to
the MS-RG stage (“g” for grazing, “pib” for the beginning of plunge-in phase,
and “piww” for when the MS star is well within the CE); according to whether
the jets were self-regulated (sr), constant (c), or if no jets were present
(nj); and according to the $\eta$ value. The models which have successful jets
are indicated by a check-mark.

Model | Stage | Jets |  $\eta_{\rm sr}$ or $\eta_{\rm c}$ | Outcome
---|---|---|---|---
g-nj-0.00 | grazing | NJ | 0.00 | $X$
g-sr-0.02 | grazing | SR | 0.02 | $X$
g-sr-0.05 | grazing | SR | 0.05 | $X$
g-sr-0.10 | grazing | SR | 0.10 | $\checkmark$
g-sr-0.30 | grazing | SR | 0.30 | $\checkmark$
g-sr-0.50 | grazing | SR | 0.50 | $\checkmark$
g-c-0.02 | grazing | C | 0.02 | $\checkmark$
g-c-0.05 | grazing | C | 0.05 | $\checkmark$
pib-nj-0.00 |  PI-beginning | NJ | 0.00 | $X$
pib-sr-0.02 |  PI-beginning | SR | 0.02 | $X$
pib-sr-0.05 |  PI-beginning | SR | 0.05 | $X$
pib-sr-0.10 |  PI-beginning | SR | 0.10 | $X$
pib-sr-0.30 |  PI-beginning | SR | 0.30 | $X$
pib-sr-0.50 |  PI-beginning | SR | 0.50 | $\checkmark$
pib-c-0.02 |  PI-beginning | C | 0.02 | $X$
pib-c-0.05 |  PI-beginning | C | 0.05 | $X$
pib-c-0.10 |  PI-beginning | C | 0.10 | $X$
pib-c-0.30 |  PI-beginning | C | 0.30 | $\checkmark$
pib-c-0.50 |  PI-beginning | C | 0.50 | $\checkmark$
piww-nj-0.00 |  PI-well within | NJ | 0.00 | $X$
piww-sr-0.05 |  PI-well within | SR | 0.05 | $X$
piww-sr-0.10 |  PI-well within | SR | 0.10 | $X$
piww-sr-0.30 |  PI-well within | SR | 0.30 | $X$
piww-sr-0.50 |  PI-well within | SR | 0.50 | $X$
piww-sr-0.70 |  PI-well within | SR | 0.70 | $X$
piww-c-0.05 |  PI-well within | C | 0.05 | $X$
piww-c-0.10 |  PI-well within | C | 0.10 | $X$
piww-c-0.30 |  PI-well within | C | 0.03 | $X$
piww-c-0.50 |  PI-well within | C | 0.50 | $X$
piww-c-0.70 |  PI-well within | C | 0.70 | $X$
Table 1: Stage and jet characteristics of each model. The models are labels
according to: a) the MS-RG stage (g for grazing, pib for when the plunge-in is
beginning, and piww for when the MS star is well within the CE), b) if the
jets were self-regulated (sr), or constant (c), or if no jets were present
(nj), and c) the efficiency value. The used acronyms are PI for plunge-in, SR
for self-regulated, C for constant. The models which have successful jets are
indicated by a check-mark, choked jets (or no jets) are indicated by $X$.

##  3 Results

In this section we show the stellar envelope morphology and its evolution as a
result of a MS star which launches two jets perpendicular to the orbital plane
either when it grazes or when it plunges-in through the envelope of a RG. We
also discuss the differences between having a pair of self-regulated jets vs
constantly powered jets.

###  3.1 Grazing stage

We first describe the expected outcomes from when the envelope of the RG star
is grazed by a MS star. The MS star can launch jets with either enough ram
pressure to overcome the stellar envelope material which surrounds the MS star
and to break out of the CE (which we term as “successful"), or jets which do
not overcome the ram pressure of the envelope material and are not able to
break out of the CE (“choked"). In column 5 of Table 1 we list whether the
jets from each model are successful or choked at the end of each simulation.

We find that in the grazing stage successful jets are easily launched (they
require $L_{j}\sim 10^{37}$ erg s-1), and no clear disk forms around the MS
star. Figure 2 shows the results for model with $\eta=0.10$ (g-sr-0.10). The
meridional (upper panels) and equatorial (lower panels) planes with the
density-map plots and the velocity field at three different times are shown.
The left panels show the initial condition ($t=t_{0}+0$ s, with $t_{0}=0$
days) when the RG is grazed by the MS star and the jets are about to be
launched. The envelope of the RG has density values which range from $\sim
10^{-8}$ g cm-3 in the outskirts of the envelope to $\sim 10^{-5}$ g cm-3 near
the core. The ambient medium around the RG is static, however since the
reference system is placed in the MS, thus, the grazing motion of the MS
around the RG star produces the initial $v_{y}\neq 0$ clearly seen in the
equatorial plane. The middle panels show a snapshot ($t=t_{0}+5\times 10^{5}$
s) when the jets are successful. Each jet propagates through a low-density
channel formed within the bulge (centered in the MS star). The bulge has
densities of order $\sim 10^{-8}-10^{-6}$ g cm-3. The jets have density and
velocity values of order $\sim 10^{-9}$ g cm-3 and $10^{7}$ cm s-1
respectively. A wind (with $\sim 10^{-8}$ g cm-3 and $v\sim 10^{6}$ cm s-1) is
produced. Meanwhile, a fraction of the material that is being accreted by the
gravitational pull of the MS star circulates in the equatorial plane around
the MS star (with $\sim 10^{7}$ cm s-1).

The panels on the right show the configuration at the end of the grazing stage
simulation ($t=t_{0}+10^{6}$ s). The jets are successful and the bulge is
still present. The density of the jets and their velocity are practically
unchanged as the jets are now further away from the MS star. The bulge
basically maintains the same properties (density, equatorial morphology flow,
and wind) and has now expanded. By the end of the integration time the initial
orbital separation ($a_{\rm orb}=R_{\rm RG}=83~{}\mbox{$R_{\odot}$}$)
diminishes to $a_{\rm orb}=0.99~{}R_{\rm RG}=82.4~{}\mbox{$R_{\odot}$}$. The
self-regulated jets with $\eta>0.10$ had enough ram pressure for the jets to
also be successful and follow the same global evolution.

Figure 2: Density map plots showing the evolution of a successful self-
regulated (SR) pair of jets with $\eta$ = 0.10 (model g-sr-0.10), when the MS
star is grazing the envelope of the RG. The meridional (upper) and equatorial
(lower) planes are shown at characteristic times ($t=t_{0}+0s$, $t_{0}+5\times
10^{5}s$, $t_{0}+10^{6}s$, with $t_{0}=0$ days). The initial orbital
separation at $t_{0}$ is $a_{\rm orb,init}=R_{\rm RG}$. The axis are in units
of $10^{12}$ cm. The velocity field is indicated by the arrows (white:
v=$10^{5}$ cm s-1, grey: v=$5\times 10^{5}$ cm s-1, black: v=$10^{7}$ cm s-1).
The magnitude of the component of the velocity in the equatorial plane is
indicated by the isocontours (black: $v_{\rm eq}=6.0\times 10^{6}$ cm s-1,
grey: $v_{\rm eq}=8.6\times 10^{6}$ cm s-1, and white: $v_{\rm eq}=1.7\times
10^{7}$ cm s-1). The left panels show the resolution at $t=t_{0}$ taken from
S19 (with $\delta X\sim$7$\times 10^{10}$ cm). As soon as $t>t_{0}$, the AMR
adapts to finer resolutions ($\delta X\sim$3.44$\times 10^{9}$ cm). The
envelope was stable without the presence of the MS and increased its size by
<5% of its initial value during the integration time.

In Figure 3 we show the evolution of the MS-RG system during the grazing stage
when no jets (NJ) are launched (model g-nj-0.00, two left panels) and when the
pair of successful jets are constantly (C) powered (with $\eta$ = 0.02, model
g-c-0.02, two right panels). The meridional and equatorial planes with the
density-map plots and the velocity field at two different times are shown
($t_{0}+5\times 10^{5}$ s and $t_{0}+10^{6}s$, the initial setup is not shown
as it is the same as for the self-regulated case). The two left panels show
how in the NJ model the bulge forms without the presence of a jet and is
basically the same as that of the SR model shown in Figure 2 (as well as the
accreted material which circulates around the MS star in the equatorial
plane). A low-density channel is still formed within the bulge, but now with
material falling towards the MS star. The two right panels show the evolution
of a pair of successful and constantly powered jets (model g-c-0.02). The
bulge morphology, the jet density and velocities are similar to those of the
SR models. The jets initially move nearly vertically (see the zoomed inset in
the upper-right panel). As the jets interact with the bulge, they are
deflected and flow through the low-density channel (see the velocity field).
Constantly-powered jets with $\eta>0.02$ present the same global morphology,
evolution, and are also successful. The equatorial flow morphology is not
affected by the presence or nature of the jets (SR or C) given that the jets
are perpendicular to the XY axis.

Figure 3: Density map plots showing the evolution during the grazing stage
when no jets (NJ) are launched (model g-nj-0.00, two leftmost panels) and when
the pair of successful jets are constantly powered (model g-c-0.02, two
rightmost panels). In each case the configuration at the characteristic times
$t_{0}+5\times 10^{5}$ s and $t_{0}+10^{6}$ s (with $t_{0}=0$ days) are shown.
The initial setup ($t=t_{0}+0s$) is not shown as it is the same as for the
self-regulated case (left panel in Figure 2). The axis units, $a_{\rm
orb,init}$, velocity field magnitude, and isocontour values are the same as in
Figure 2. The inset in the upper-right panel shows a zoom of the inner region
where the jets are being launched.

The mass-accretion rate and specific angular momentum for the successful jet
models and the choked (or NJ) models, during the grazing stage ($t=t_{0}$),
are shown in Figure 4 and Figure 5, respectively. Specifically, we show the
mass accretion and specific angular momentum at $R_{\rm{in}}$
($\dot{M}_{\rm{in}}$, $J_{\rm sp,in}$, respectively) for the successful SR
model g-sr-0.10 (termed as SR high $\eta$), the successful C model g-c-0.02 (C
high $\eta$), and the NJ (model g-nj-0.00), i.e. SR and C models with low
$\eta$ which choke (NJ or SR-C low $\eta$). The $\dot{M}_{\rm{in}}$ of the
successful self-regulated jets tends to be marginally lower than that of the
successful constantly-powered jets. The choked-jet models tend to have
slightly higher accretion rates. The mass accretion rate is nearly identical
in all the cases (this is, either when no jets were launched, when they were
choked, or when they were successful). Independently of the nature of the jets
(choked or not), the mass accretion at the inner boundary was within
$\dot{M}_{\rm{in}}\approx 2-6\times 10^{23}$ g s-1. Thus, the self-regulated
jets during the grazing stage had a $L_{j}\sim 10^{36}-10^{37}$ erg s-1
(depending on the $\eta$ value), which is about an order of magnitude below
the constantly-powered-jet models luminosity. As for the accretion rate, in
the choked jet models the $J_{\rm sp,in}$ tends to be moderately higher than
that of the successful jets. Still, there is no clear difference in the
$J_{\rm sp,in}$ between the successful self-regulated jets, the constantly
powered successful jets, and the choked jet models. For all cases $J_{\rm
sp,in}\approx 0.4-1.4\times 10^{18}$ cm2 s-1.

Figure 4: Mass accretion rates at $R_{\rm in}$ during: the grazing stage, $a)\
t=t_{0}$; the beginning of the plunge-in stage, $b)\ t=t_{pib}$; and when the
MS star is well within the CE, $c)\ t=t_{piww}$. Self-regulated jet models
with high $\eta$ are indicated in blue, constantly powered jets with high
$\eta$ in red, and low $\eta$ choked jets in black. The horizontal dashed
lines only separate cases a), b), and c) and have no physical meaning. Figure
5: Specific angular momentum at $R_{\rm in}$ for the models indicated in
Figure 4.

###  3.2 Plunge-in stage

We next discuss the global evolution of the system as the MS star plunges-in
through the envelope of the RG. The jets are launched once the MS is within
the CE, i.e., either during the pib or piww stages. Table 1 lists whether the
jet models at each stage are successful or choked at the end of the
integration time.

In the large scale global simulation of S19 which is used as our initial
configuration, the plunging-in of the MS star heats the envelope at the MS
star surroundings. As a result, the envelope is inflated (and remains mostly
bound), and an equatorial outflow is produced. In addition, envelope mass
accumulates behind the path of the MS and concentrates around the equatorial
plane in a bulge-like structure. Later, when the MS star further plunges-in,
both the MS star and the bulge are engulfed inside the envelope (see left-most
panels of figures 6 and 7). Because initially the envelope in the companion’s
vicinity concentrates around the equatorial plane the jets will have less
layers of envelope to interact with in the vertical direction, which will ease
the launching of successful jets, while in the piww stage when the bulge is
immersed inside the envelope there is a greater chance that the jets would be
choked.

In the pib stage successful jets require higher $\eta$ powering (compared to
the grazing stage). Figure 6 shows the meridional and equatorial density maps
and velocity field evolution (upper and lower panels, respectively) for the SR
model with $\eta=0.50$ during the pib stage (model pib-sr-0.50). The orbital
separation in this stage changes from $a_{\rm orb}=0.76~{}R_{\rm RG}$ to
$a_{\rm orb}=0.49~{}R_{\rm RG}$. The three leftmost panels show the initial
condition and two subsequent times ($5\times 10^{5}$ s and $10^{6}$ s)
exhibiting how the jets are able to drill through the RG. For comparison, a
snapshot of the model where no jets are launched is also included (pib-nj-0.00
at $5\times 10^{5}$ s, rightmost panel). The morphology of the bulge and the
region accreting onto the secondary star are very similar in the cases
with/without jets. The bulge, composed once more by a wind and equatorial
material (with densities $\sim 10^{-6}$ g cm-3), is present, and no clear disk
forms around the MS star. The jets have density and velocity values of order
$\sim 10^{-8}$ g cm-3 and $10^{7}$ cm s-1 respectively. Since the material
that surrounds the MS star is denser than in the grazing stage, the jets must
be more powerful in order to be successful. For the self-regulated jets to be
successful the $\eta$ efficiency must be much larger ($\eta\geq 0.50$) than
that in the grazing stage. Thus, the jets now require a luminosity $\sim
10^{39}$ erg s-1 in order to be successful. Constantly powered jet models with
$\eta>0.30$ value are also successful, and present a similar global morphology
and evolution.

In the piww stage all jets, independently of their nature or $\eta$
efficiency, are choked. Figure 7 shows the meridional and equatorial density
maps and velocity field evolution (upper and lower panels, respectively) for
the SR model with $\eta=0.70$ during the piww stage (model piww-sr-0.70). The
orbital separation in the piww stage changes from $a_{\rm orb}=0.48~{}R_{\rm
RG}$ to $a_{\rm orb}=0.47~{}R_{\rm RG}$. The three leftmost panels show the
initial condition and two subsequent times when the MS is well within the RG.
Also, a snapshot of the model where no jets are launched is also included
(piww-nj-0.00, rightmost panel). The high density which surrounds the MS star
is such that the ram pressure of the jet is not able to drill through the core
of the RG, and thus the jets are always choked (independently if they are SR
or C powered). The jets manage to initially expand outside the inner boundary
but are quickly choked ($\sim 10^{5}s$) inside the bulge (which is slightly
modified by the presence of the choked jet). The equatorial morphology flow is
not affected by the presence of jets and no disk forms around the MS star.

Figure 6: Density map plots showing the evolution of a successful SR jets with
$\eta$ = 0.50 (model pib-sr-0.50) for when the MS star is starting the plunge-
in phase through the envelope of the RG (where $t_{pib}=35$ days, and $a_{\rm
orb,init}=0.76R_{\rm RG}$). The meridional (upper) and equatorial (lower)
planes are shown at characteristic times ($t=t_{pib}+0s$, $t_{pib}+5\times
10^{5}s$, $t_{pib}+10^{6}s$). A snapshot of model pib-nj-0.00 is also included
(upper and lower right panels). The axis units, velocity field magnitude, and
velocity isocontours are the same as Figure 2. Figure 7: Same as Figure 7 for
when the MS star is well within the RG ($t_{piww}=52$ days, nd $a_{\rm
orb,init}=0.48R_{\rm RG}$).

The $\dot{M}_{\rm{in}}$ and $J_{\rm sp,in}$ for the successful SR or C jet
models and the NJ or choked models, during the pib and piww stage
($t=t_{\rm{pib}}$, and $t=t_{\rm{piww}}$), are shown in Figure 4 and Figure 5,
respectively. In the pib stage the we show model pib-sr-0.50 as the successful
SR model, model pib-c-0.30 as the successful C model, and pib-nj-0.00 for the
choked model. In the piww stage the we show models piww-sr-0.70, piww-c-0.70,
and piww-nj-0.00. During the pib the accretion rates are within
$\dot{M}_{\rm{in}}\approx 0.7-4.0\times 10^{24}$ g s-1, which is about an
order of magnitude larger than in the grazing stage. Meanwhile, at the piww
stage, where all the jets were choked, the accretion rate was $\sim$ two
orders of magnitude larger than in the grazing stage, this is,
$\dot{M}_{\rm{in}}\approx 0.7-3.0\times 10^{25}$ g s-1. The $J_{\rm sp,in}$
during the pib stage is $J_{\rm sp,in}\approx 0.4-1.7\sim 10^{18}$ cm2 s-1,
and at the piww stage is $J_{\rm sp,in}\approx 1.0-1.8\sim 10^{18}$ cm2 s-1.

##  4 Discussion

We have studied the dynamics and propagation of jets launched from a MS star
moving through the envelope of a RG star. We followed, using 3D HD
simulations, a set of jet models (either self-regulated or constantly powered,
and with different kinetic luminosities) in three phases: when the MS is
grazing the RG, when the MS star has just started to plunge-in through the CE,
and when the MS star is well within the CE. Our numerical simulations show
that jets can be choked or successful depending on the jet efficiency and on
the evolutionary phase. Grazing jets and jets launched at the beginning of the
plunge-in stage may be successful. Meanwhile, jets that are launched deep
inside the CE are choked. As shown in Section 3, it becomes increasingly
difficult for the jets to break out through the accreting material as the MS
star moves inwards into the CE. Once the jets are launched, independently of
whether they are self-regulated or constantly powered, each jet will be
successful if its ram pressure ($P_{\rm jet}$) is larger than the ram pressure
of the accreting ambient material ($P_{\rm amb}$). The jet and the accreting
material ram pressures are given by (see Moreno Méndez, López-Cámara, & De
Colle 2017 for further details)

$\displaystyle P_{\rm jet}\simeq\rho_{\rm
amb}v_{\star}^{2}\frac{r_{B}^{2}}{r^{2}}\frac{2\eta
v_{j}}{\theta_{j}^{2}v_{\star}}\;,\qquad P_{\rm amb}\simeq\rho_{\rm
amb}v_{\star}^{2}\frac{r_{B}^{2}}{r^{2}}\;,$ (1)

where $\rho_{\rm amb}$ is the density of the CE, $v_{\star}=\sqrt{GM(a)/a}$ is
the velocity of the MS star within the RG (specifically at the orbital
separation $a$, and taken as keplerian), $r_{B}=2GM_{\rm MS}/v_{\star}^{2}$ is
the Bondi radius, $v_{j}$ the velocity of the jet, and $\theta_{j}$ the
opening angle of the jet. From the ram pressure of the jets we can see that if
they were faster, or if they were more collimated then smaller $\eta$
efficiencies would be required for the jets to be successful. More collimated
jets have a much larger ram pressure (since $P_{\rm
jet}\propto\theta_{j}^{-2}$), which may overcome the accreting material much
more easily. In our case, if the opening angle of the jets was
$\theta_{j}\approx 10^{\circ}$ (instead of the $\theta_{j}\approx 30^{\circ}$
used in our simulations) then the $\eta$ values for the jets to be successful
would be $\sim$10 times smaller.

Comparing the ram pressures, the condition for the jets to be successful
$P_{\rm jet}>P_{\rm amb}$ implies that

$a\gtrsim a_{0}=3\times 10^{12}{\rm
cm}\left(\frac{M(a)}{1M_{\odot}}\right)\left(\frac{\theta_{j}}{30^{\circ}}\right)^{4}\left(\frac{\eta}{0.1}\right)^{-2}\left(\frac{v_{j}}{500\;{\rm
km}\;{\rm s}^{-1}}\right)^{-2}.$ (2)

Equation 2 shows that outside a critical orbital separation $a_{0}$, jets are
successful, while at smaller orbital separations the jets will be choked by
the accreting material. The larger the efficiency, hence, smaller the orbital
separation (and deeper within the CE) where the jets are successful. We also
see that the value of $a_{0}$ strongly depends on the jet velocity, with
larger velocity the jet ram pressure increases, thus making it easier for a
jet to break out from deeper within the CE. The opening angle also noticeably
modifies the critical orbital separation in which the jets would be successful
(as $a_{0}\propto\theta_{j}^{4}$). Hence, jets with small efficiency could
produce successful jets if their opening angle was much smaller, at least in
the intermediate “plunge-in” phase. In order for the jets to be successful,
their ram pressure ought to be larger than that of the surrounding stellar
envelope (both computed in the shock frame). This is true both in accreting
systems and in static (not-accreting) flows, in which the ram pressure of the
ambient medium corresponds to its inertia. The thermal pressure is not large
enough to slow the jet down as the jet shock head is highly supersonic and in
strong shocks the thermal pressure of the ambient medium is negligible.

As the jet velocity is typically close to the escape velocity,
$v_{j}\propto(M/R)^{1/2}$ (where $M$ and $R$ are the mass and radius of the
object which is launching the jet), then equation 2 leads to $a_{0}\propto
R/M$. Since for MS stars, $R\propto M^{0.8}$ (e.g. Popper, 1980), then $a_{0}$
depends weakly on the mass of the main sequence star ($a_{0}\propto
M^{-0.2}$). Thus, the results obtained in this paper, in particular where the
jet is successful or choked, can be generalized to low mass stars ($M\lesssim
1$ $M_{\odot}$). A very different outcome is expected in the case of neutron
stars and black holes in which the radius is several orders of magnitude
smaller than that for MS stars, thus, $a_{0}$ is much smaller and jets are
energetic enough to propagate successfully at any orbital separation (Moreno
Méndez, López-Cámara, & De Colle, 2017). The accretion rate in all stages was
computed by: i) considering the enclosed mass in the orbital separation using
the initial density profile (and assuming spherical symmetry to calculate
$v_{kep}$), and ii) using the density from the original profile close to the
location of the MS star at each stage. The accretion rate found in the grazing
stage is $\sim 3\times 10^{23}$ g s-1, in the pib $\sim 2\times 10^{24}$ g
s-1, and in the piww $\sim 10^{25}$ g s-1. The accretion rates we find in all
stages are approximately 1-10% of the correspondent BHL accretion rate. The
latter is consistent with the mass accretion rate onto a neutron star (or a
black hole) immersed within the CE of a massive star (MacLeod & Ramirez-Ruiz,
2015; Moreno Méndez, López-Cámara, & De Colle, 2017; Gilkis, Soker, & Kashi,
2019; López-Cámara, De Colle, & Moreno Méndez, 2019; López-Cámara, Moreno
Méndez, & De Colle, 2020). Thence, the successful jets during the grazing and
beginning of the plunge-in stage require luminosities of order $L_{j}\gtrsim
10^{37}$ erg s-1 and $L_{j}\gtrsim 10^{39}$ erg s-1, respectively. When the MS
star is well within the CE, according to equation 2, the jets are choked even
for $10^{40}$ erg s-1. Also, the successful self-regulated jet models need
less power than the constantly powered models (by a factor of $\sim$2-3). For
the jets to be successful they need to clear the region close to the injection
boundary. This is done more efficiently if the jets self-regulate and their
luminosity adjusts to the changing conditions of the accreting environment.

The Eddington luminosity is $L_{\rm{Edd}}\sim 10^{38}(M/M_{\odot})$ erg s-1,
while the release of the binding energy of the accreting material onto the
stellar surface is five to a hundred times larger. Thus, we have super-
Eddington luminosity powering the jets. We notice that to properly determine
the fate of the material crossing the inner boundary requires simulations
which include radiation transport (due to the super-Eddington accretion) and
in which the inner boundary is set at the stellar surface of the MS star
(which is a factor of $\sim$ 3 smaller than our inner boundary).

We have assumed that the jets from each model are emitted at a stage of the MS
and RG system evolution (either during the grazing-stage, or at the
commencement or termination of the plunge-in stage). The delay in the
launching of the jets may be due to the fact that an accretion disk around the
MS star may not have formed, or that the jet launching mechanism comes into
play. The disk can be formed through two main channels, either by the
accretion of the CE material, or if the disk was pre-existent (i.e., formed
before the MS star was engulfed by the RG and survived the grazing and plunge-
in phase). The critical specific angular momentum to form a centrifugally
supported disk at the inner boundary is $J_{\rm kep,in}=R_{\rm in}v_{\rm
kep,in}$, with $v_{\rm kep,in}$ the Keplerian velocity at the inner boundary
($\approx 2\times 10^{7}$ cm s-1). Thus, $J_{\rm kep,in}=R_{\rm in}v_{\rm
kep,in}\approx 2\times 10^{18}$ cm2 s-1. Since the velocity in the orbital
plane at $R_{\rm in}$ is $\approx 10^{7}$ cm s-1, then the specific angular
momentum at the inner boundary is $J_{\rm sp,in}\approx 10^{18}$ cm2 s-1. The
latter corresponds to $\sim$0.5$J_{\rm kep,in}$, thus, our calculations show
that it is unlikely that all accreted material will be accreted through a disk
but a fraction of it possibly will. Our simulations lack magnetic fields and
other angular momentum transfer mechanisms (for example tidal synchronization,
shear, viscosity, etc, see Armitage (2022) for further details), thus our
angular momentum values should be considered as un upper limit.

The MS star is expected to enter the common envelope with a disk from a pre-
phase of RLOF. Such disks are observed in several common envelope simulations,
e.g., S19, Staff et al. (2016), Reichardt et al. (2019). The question is
whether such disks extend to the smaller scale around the companion and
whether accretion can take place at such disks. For this purpose, simulations
that include accretion modeling and that resolve the companion need to be
carried-out. In our simulations, using a sub-grid accretion model with which
we are close to resolving the MS star we find that the critical angular
momentum to form a centrifugally supported accretion disk at the surface of a
$0.3M_{\odot}$ star ($R_{\rm surf}\sim 3\times 10^{10}$ cm) is $J_{\rm
kep,surf}\simeq 10^{18}$ cm2 s-1. If the angular momentum is conserved between
$R_{\rm in}$ and $R_{\rm surf}$, then our simulations indicate that a small
centrifugally supported disk may form close to the surface of the MS star
(since the specific angular momentum may be larger than $J_{\rm kep,surf}$).
In order to properly address the disk creation around the MS simulations with
a smaller $R_{\rm in}$ than that in our study, simulations which include
radiation transport and magnetic fields (which properly transport the angular
momentum, Balbus & Hawley, 1991) are required. Given that binary systems that
undergo a CE phase have to previously experience a Roche-lobe overflow (RLOF)
mass-transfer phase, thus, it is likely that a accretion disk may have formed
around the MS star before it was engulfed by the CE. The latter has been
observed in CE simulations that start with a companion far from the RG ($a\geq
2R_{\rm RG}$) and in which a disk is formed around the MS star (in particular
when the pre-RLOF stage is simulated Shiber & Iaconi, 2021). It is also
possible that the sub-Keplerian inflows will form an accretion belt around the
accreting object (instead of an accretion disk) from which jets may be
launched (Schreier & Soker, 2016), or a polar outflow may be produced without
the presence of an accretion disk due to an inflow with an equatorial-to-polar
density gradient (Lery et al., 2002; Aguayo-Ortiz, Tejeda, & Hernandez, 2019;
Aguayo-Ortiz, Sarbach, & Tejeda, 2021).

It is unclear what would be the long-term evolution of the system in the case
of choked jets. The continuous injection of energy from the star/disk system
into the jets will be deposited in the region close to the star. If the ram
pressure of the accreted material is much larger than the jet ram pressure,
the jet material flow will be reversed and it will go back to the disk,
forming a small “circuit” of ejected-accreted material (in which most likely
the accreted jet material will be a fraction of the accreted envelope
material). Detailed numerical simulations are required to compute this
precisely.

The setups used as the initial condition of our simulations (at the grazing,
pib and piww stages) taken from S19 did not consider jets. Nevertheless, the
difference in the orbital separation and in the mass ejection when considering
jets is negligible (though the geometrical envelope structure is somewhat
different). Thus, our results can be compared with the large scale models of
S19 that do include jets. Both studies agree on the fact that, in the grazing
stage, even jets with low $\eta$-efficiencies may be successful. Also the jets
are diverted outwards by the envelope of the RG, thus, the resulting
geometrical structure produced around the edge of the envelope and the MS star
is very similar. At the pib stage, the jets of the large scale simulation
successfully propagate outside, while in the local simulation they are only
successful for relatively large $\eta$ values. In the piww stage the jets are
initially choked inside the envelope in both simulations. However, in the
large scale simulation the jets are able to drill outside on a time scale of 3
days, whereas in the small scale simulations they are still choked after $\sim
10$ days.

Based on the luminosity estimates of order $L_{j}\gtrsim 10^{37-39}$ erg s-1
and integrating over jets energy deposition timescales of order
$10^{2}-10^{3}$days (as obtained in the 3D HD large scale simulations of S19),
we have total energy deposition between $E_{j}\gtrsim 10^{44-47}$ ergs.
Considering that the binding energy of the RG star is $E_{\rm bind}\simeq
10^{47-48}$ ergs, hence, most of the CE will remain bound. However, since the
binding energy of the outer part of the envelope of the RG (i.e. that for
$R>20R_{\rm\odot}$) is $\sim 10^{46}\;{\rm ergs}$, thus, jets may facilitate
partial CE-envelope removal (including other sources of energy such as
recombination, magnetic fields, enthalpy could help to fully remove the CE,
see Ivanova et al. (2013) for further details). A consequence of the partial
CE-stripping would be that less material is present during the plunge-in
stage, and the jets may be successful. If the CE of the RG is mostly
convective, which would be the case when the low mass RG has evolved from its
respective terminal main sequence stage, then the material is less bound
(Shore et al., 1994; Klencki et al., 2021). Thus, the CE would be lifted more
efficiently by the jets.

The total amount of accreted mass in our simulations is still relatively high.
On the short time scale we simulated ($10^{6}\;{\rm s}\approx 11.6$ days), the
MS star accretes $\sim 10^{-2}M_{\odot}$ while plunging-in. This is $\sim$1%
of the total mass we had initially in the grid during the correspondent stage.
If such an accretion rate persists for a typical duration of $\sim 500$ days
in which the envelope has not been removed yet, the MS star can accrete $\sim
0.5M_{\odot}$. Assuming a similar amount of accretion in a CE event while the
primary is an AGB star, it is possible that the carbon enrichment of dwarf
carbon stars, which have been recently observed in short-period binary systems
of less than a day (Roulston et al., 2021; Whitehouse et al., 2021), can take
place also during the CE.

The photospheric temperature of the CE is $T_{\rm eff}\sim 10^{3}$ K, thus,
its emission properties will likely resemble those of a RG and be observable
in the NIR-optical(see for example: Valenti, Ferraro, & Origlia, 2004; Huang
et al., 2020). As the medium is optically thick, propagation of the jets
through the CE may produce a bright quasi-thermal spectrum. Jets powered by
high mass accretion during the GEE could produce a transient event since a
non-negligible fraction of the kinetic energy of the jets may be be
transferred to radiation energy (Soker, 2016). While a full radiation HD
calculation is needed to properly compute the observational outcomes of the
jet-cocoon system breaking out of the CE, as well as those of the ejected
stellar material, bright thermal emission is expected. From the Rankine-
Hugoniot jump conditions the temperature of the post-shock material can be
estimated by $kT=3/16m_{p}v_{\rm sh}^{2}$, where $v_{\rm sh}$ is the shock
velocity. Thus, for shock velocities or order $v_{\rm sh}\sim 10^{7}-4\times
10^{7}$ cm s-1 the temperature is $\sim 10^{5}-10^{6}$ K. As the region close
to the MS star is opaque (as its optical depth is $\tau\sim\rho kR\gg 1$), the
emission will be thermal and its luminosity555The large uncertainty in the
luminosity is direct consequence in the uncertainty in the temperature in the
photosphere (which is not spatially resolved in our simulations). will be of
order $L=4\pi\sigma R^{2}T^{4}\sim(10^{39}-10^{43})\ R_{11}^{2}$ erg s-1.
Thus, the spectrum will peak in X-rays and UV at early times and evolves to
lower frequencies as a function of time. Transients with these luminosities
are visible up to distances of $\approx$ 1-100 Mpc.

The jets injection rate of simulation #6 in S19 was 0.001 M⊙ yr-1, which was
0.66% of the BHL accretion rate calculated at the RG surface $\sim$0.15 M⊙
yr${{}^{-}1}$. However, we find much lower accretion rates in the grazing (30
times smaller) and pib (5 times smaller) stages, and about the same accretion
rates in the piww stage. This would imply efficiency values of 20%, 3.5%, and
0.7% of the jets in S19, at each stage, respectively. Therefore, it might be
that the jets power in S19 is overestimated in the grazing and pib phases,
although these efficiency values are still reasonable considering the fact
that power of the jets in S19 was constant through the simulation evolution.
As the main effect of the jets on the orbital evolution and on the envelope,
unbinding starts only afterwards (t>53 days), this overestimation most likely
did not play a significant role.

On the other hand, because most of the unbinding in S19 happens afterwards
when the MS is deep inside the envelope, the question whether the jets are
successful to penetrate through the envelope or choked at this stage is
critical. In our simulation the jets are choked well within the envelope and
therefore are not expected to be able to unbind the envelope. There are two
factors which make a difference on whether the jets propagate or not in the
simulations of S19 and the ones in this paper: 1) The injection radius in our
case is 1011 cm, while in S19 is 1012 cm. Thus, our jets have to overcome a
larger amount of ambient medium with a higher ram pressure. The energy
necessary to move the extra mass, from deeper within the gravitational
potential, is not considered in simulations that launch the jets further
outside. 2) The boundary conditions responsible for the jet injection in S19
are different with respect to our paper. In both simulations the jets are
injected from a cone with an opening angle of 30∘. While in S19 the density is
redefined at each timestep, in our simulations both density and pressure are
redefined at each timestep. To better address this question a future study
that explores in detail the numerical parameters involved in the simulations
is required.

##  5 Conclusions

In this paper, we have presented small scale 3D numerical simulations of the
propagation of jets through a CE formed by a 0.88 $M_{\odot}$ red giant and a
0.3 $M_{\odot}$ main sequence star. The simulations were based on the CE
structure configuration resulting from the large scale simulations of Shiber
et al. (2019). We find that the final outcome depends on the conditions of the
environment, on the power of the jet, and whether the jet is constantly
powered or self-regulated. In the grazing stage and the commencement of the
plunge self-regulated jets need higher efficiencies to break out of the
envelope of the RG. Deep inside the CE, on the timescales simulated, jets are
choked independently of whether they are self-regulated or constantly powered.
Jets with the same characteristics as those in the large-scale simulations of
S19 (which are able to break out of the envelope of the RG), are choked in our
small-scale simulations (for at least the integration time, $\sim$11 days).
Having a smaller jet injection radius, and the self-regulation of the jets may
affect the outcome. Our simulations are a step in the direction of
investigating the interface between local and global simulations.

Once jets are launched, their success and propagation through the CE depends
on the conditions of the environment as well as on the power of the jets.
Self-regulated jets will be by their nature more efficient than the constantly
powered case. These are more luminous when the accretion rate is high, and
less luminous when the accretion rate is small. On the other hand, constantly
powered jets are always using the same energy deposition rate, regardless of
the accretion rate. Also, high luminosity emission is expected in this case
going from X-rays to UV and optical as a function of time. Deep inside the CE
the jets are choked since the ram pressure of the jets is overcome by that of
the accreting material.

Our simulations show that the mass accretion onto the MS star is 1-10% of the
Bondi Hoyle Littleton rate ($\sim 10^{-3}-10^{-1}$ M⊙ yr-1), and that a disk
around the secondary MS star is probably not formed. Still, jets could be
launched in the case of a pre-existing disk or via MHD effects. It is unclear
if the disk can survive when it is engulfed into the CE and needs to be
further studied. Future surveys (including the Vera C. Rubin Observatory
666www.lsst.org) will help to decipher the CE phase. The outlook for the study
of the CE phase is promising as it is likely that many of the observed HEAP
events will allow to place constraints on the models.

We must note that we only follow $\sim$11 days in each of the three stages,
thus we do not follow the totality of the interaction of jets within the CE.
Future large scale simulations with high-resolution, which are currently at
the edge of our computational capabilities, will permit to fully study the
evolution of jets, as well as the disk creation, when immersed in the CE of a
massive star.

## Acknowledgments

We thank the referee for the helpful comments and improvement of the
manuscript. D.L.C. is supported by Cátedras CONACyT at the Instituto de
Astronomía, UNAM. We acknowledge the support from the Miztli-UNAM
supercomputer (projects LANCAD-UNAM-DGTIC-321 and LANCAD-UNAM-DGTIC-281) for
the assigned computational time in which the simulations were performed.
D.L.C. and F.D.C. thank the UNAM-PAPIIT grant IG100820. S.S. thanks National
Science Foundation Award 1814967 for its support. R.I. is grateful for the
financial support provided by the Postodoctoral Research Fellowship of the
Japan Society for the Promotion of Science (JSPS P18753). Many of the images
in this study were produced using VisIt. VisIt is supported by the Department
of Energy with funding from the Advanced Simulation and Computing Program, the
Scientific Discovery through Advanced Computing Program, and the Exascale
Computing Project.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, PhRvL, 116, 241103
  * Aguayo-Ortiz, Sarbach, & Tejeda (2021) Aguayo-Ortiz A., Sarbach O., Tejeda E., 2021, PhRvD, 103, 023003. doi:10.1103/PhysRevD.103.023003
  * Aguayo-Ortiz, Tejeda, & Hernandez (2019) Aguayo-Ortiz A., Tejeda E., Hernandez X., 2019, MNRAS, 490, 5078. doi:10.1093/mnras/stz2989
  * Armitage (2022) Armitage P. J., 2022, arXiv, arXiv:2201.07262
  * Armitage & Livio (2000) Armitage P. J., Livio M., 2000, ApJ, 532, 540. doi:10.1086/308548
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214. doi:10.1086/170270
  * Bear, Grichener, & Soker (2017) Bear E., Grichener A., Soker N., 2017, MNRAS, 472, 1770. doi:10.1093/mnras/stx2125
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883. doi:10.1093/mnras/199.4.883
  * Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273. doi:10.1093/mnras/104.5.273
  * Brown et al. (2007) Brown G. E., Lee C.-H., Moreno Méndez E., 2007, ApJL, 671, L41
  * Castelletti et al. (2006) Castelletti G., Dubner G., Golap K., Goss W. M., 2006, A&A, 459, 535. doi:10.1051/0004-6361:20065061
  * Chamandy et al. (2020) Chamandy L., Blackman E. G., Frank A., Carroll-Nellenback J., Tu Y., 2020, MNRAS, 495, 4028. doi:10.1093/mnras/staa1273
  * Chamandy et al. (2019) Chamandy L., Blackman E. G., Frank A., Carroll-Nellenback J., Zou Y., Tu Y., 2019, MNRAS, 490, 3727. doi:10.1093/mnras/stz2813
  * Chamandy et al. (2018) Chamandy L., Frank A., Blackman E. G., Carroll-Nellenback J., Liu B., Tu Y., Nordhaus J., et al., 2018, MNRAS, 480, 1898. doi:10.1093/mnras/sty1950
  * Chevalier (2012) Chevalier R. A., 2012, ApJL, 752, L2
  * Clayton et al. (2017) Clayton M., Podsiadlowski P., Ivanova N., Justham S., 2017, MNRAS, 470, 1788. doi:10.1093/mnras/stx1290
  * De et al. (2020) De S., MacLeod M., Everson R. W., Antoni A., Mandel I., Ramirez-Ruiz E., 2020, ApJ, 897, 130. doi:10.3847/1538-4357/ab9ac6
  * De Colle et al. (2012) De Colle F., Granot J., López-Cámara D., Ramirez-Ruiz E., 2012a, ApJ, 746, 122
  * De Marco & Izzard (2017) De Marco O., Izzard R. G., 2017, PASA, 34, e001. doi:10.1017/pasa.2016.52
  * Everson et al. (2020) Everson R. W., MacLeod M., De S., Macias P., Ramirez-Ruiz E., 2020, ApJ, 899, 77. doi:10.3847/1538-4357/aba75c
  * Gilkis & Soker (2014) Gilkis A., Soker N., 2014, MNRAS, 439, 4011. doi:10.1093/mnras/stu257
  * Gilkis, Soker, & Kashi (2019) Gilkis A., Soker N., Kashi A., 2019, MNRAS, 482, 4233. doi:10.1093/mnras/sty3008
  * Gilkis, Soker, & Papish (2016) Gilkis A., Soker N., Papish O., 2016, ApJ, 826, 178. doi:10.3847/0004-637X/826/2/178
  * Glanz & Perets (2018) Glanz H., Perets H. B., 2018, MNRAS, 478, L12. doi:10.1093/mnrasl/sly065
  * Glanz & Perets (2021a) Glanz H., Perets H. B., 2021, MNRAS, 500, 1921. doi:10.1093/mnras/staa3242
  * Glanz & Perets (2021b) Glanz H., Perets H. B., 2021, MNRAS.tmp. doi:10.1093/mnras/stab2291
  * González-Casanova et al. (2014) González-Casanova, D. F., De Colle, F., Ramirez-Ruiz, E., et al. 2014, ApJ, 781, L26. doi:10.1088/2041-8205/781/2/L26
  * Grichener & Soker (2017) Grichener A., Soker N., 2017, MNRAS, 468, 1226. doi:10.1093/mnras/stx534
  * Han et al. (1995) Han Z., Eggleton P. P., Podsiadlowski P., Tout C. A., 1995, MNRAS, 277, 1443. doi:10.1093/mnras/277.4.1443
  * Hoyle & Lyttleton (1939) Hoyle F., Lyttleton R. A., 1939, PCPS, 35, 405. doi:10.1017/S0305004100021150
  * Huang et al. (2020) Huang C. D., Riess A. G., Yuan W., Macri L. M., Zakamska N. L., Casertano S., Whitelock P. A., et al., 2020, ApJ, 889, 5. doi:10.3847/1538-4357/ab5dbd
  * Iaconi et al. (2019) Iaconi R., Maeda K., De Marco O., Nozawa T., Reichardt T., 2019, MNRAS, 489, 3334. doi:10.1093/mnras/stz2312
  * Iaconi et al. (2020) Iaconi R., Maeda K., Nozawa T., De Marco O., Reichardt T., 2020, MNRAS, 497, 3166. doi:10.1093/mnras/staa2169
  * Iaconi et al. (2017) Iaconi R., Reichardt T., Staff J., De Marco O., Passy J.-C., Price D., Wurster J., et al., 2017, MNRAS, 464, 4028. doi:10.1093/mnras/stw2377
  * Iben & Livio (1993) Iben I., Livio M., 1993, PASP, 105, 1373. doi:10.1086/133321
  * Ivanova et al. (2013) Ivanova N., Justham S., Chen X., De Marco O., Fryer C. L., Gaburov E., Ge H., et al., 2013, A&ARv, 21, 59. doi:10.1007/s00159-013-0059-2
  * Jones (2020) Jones D., 2020, rfma.book, 123. doi:10.1007/978-3-030-38509-5_5
  * Kamiński et al. (2021) Kamiński T., Steffen W., Bujarrabal V., Tylenda R., Menten K. M., Hajduk M., 2021, A&A, 646, A1. doi:10.1051/0004-6361/202039634
  * Klencki et al. (2021) Klencki J., Nelemans G., Istrate A. G., Chruslinska M., 2021, A&A, 645, A54. doi:10.1051/0004-6361/202038707
  * Khokhlov et al. (1999) Khokhlov A. M., Höflich P. A., Oran E. S., Wheeler J. C., Wang L., Chtchelkanova A. Y., 1999, ApJL, 524, L107. doi:10.1086/312305
  * Kuruwita, Staff, & De Marco (2016) Kuruwita R. L., Staff J., De Marco O., 2016, MNRAS, 461, 486. doi:10.1093/mnras/stw1414
  * Lery et al. (2002) Lery T., Henriksen R. N., Fiege J. D., Ray T. P., Frank A., Bacciotti F., 2002, A&A, 387, 187. doi:10.1051/0004-6361:20020369
  * López-Cámara, De Colle, & Moreno Méndez (2019) López-Cámara D., De Colle F., Moreno Méndez E., 2019, MNRAS, 482, 3646. doi:10.1093/mnras/sty2959
  * López-Cámara, Moreno Méndez, & De Colle (2020) López-Cámara D., Moreno Méndez E., De Colle F., 2020, MNRAS, 497, 2057. doi:10.1093/mnras/staa1983
  * MacLeod et al. (2017) MacLeod M., Antoni A., Murguia-Berthier A., Macias P., Ramirez-Ruiz E., 2017, ApJ, 838, 56. doi:10.3847/1538-4357/aa6117
  * MacLeod & Ramirez-Ruiz (2015) MacLeod M., Ramirez-Ruiz E., 2015, ApJ, 803, 41. doi:10.1088/0004-637X/803/1/41
  * Moreno Méndez et al. (2011) Moreno Méndez E., Brown G. E., Lee C.-H., Walter F. M., 2011, ApJ, 727, 29. doi:10.1088/0004-637X/727/1/29
  * Moreno Méndez, López-Cámara, & De Colle (2017) Moreno Méndez E., López-Cámara D., De Colle F., 2017, MNRAS, 470, 2929. doi:10.1093/mnras/stx1385
  * Ohlmann et al. (2016) Ohlmann S. T., Röpke F. K., Pakmor R., Springel V., 2016, ApJL, 816, L9. doi:10.3847/2041-8205/816/1/L9
  * Papish & Soker (2011) Papish O., Soker N., 2011, MNRAS, 416, 1697. doi:10.1111/j.1365-2966.2011.18671.x
  * Papish & Soker (2014) Papish O., Soker N., 2014, MNRAS, 438, 1027. doi:10.1093/mnras/stt2199
  * Passy et al. (2012) Passy J.-C., De Marco O., Fryer C. L., Herwig F., Diehl S., Oishi J. S., Mac Low M.-M., et al., 2012, ApJ, 744, 52. doi:10.1088/0004-637X/744/1/52
  * Paczynski (1976) Paczynski B., 1976, IAUS, 73, 75
  * Piran et al. (2017) Piran T., Nakar E., Mazzali P., Pian E., 2017, arXiv, arXiv:1704.08298
  * Popper (1980) Popper D. M., 1980, ARA&A, 18, 115. doi:10.1146/annurev.aa.18.090180.000555
  * Postnov & Yungelson (2014) Postnov K. A., Yungelson L. R., 2014, LRR, 17, 3
  * Prust & Chang (2019) Prust L. J., Chang P., 2019, MNRAS, 486, 5809. doi:10.1093/mnras/stz1219
  * Reichardt et al. (2020) Reichardt T. A., De Marco O., Iaconi R., Chamandy L., Price D. J., 2020, MNRAS, 494, 5333. doi:10.1093/mnras/staa937
  * Reichardt et al. (2019) Reichardt T. A., De Marco O., Iaconi R., Tout C. A., Price D. J., 2019, MNRAS, 484, 631. doi:10.1093/mnras/sty3485
  * Ricker & Taam (2012) Ricker P. M., Taam R. E., 2012, ApJ, 746, 74. doi:10.1088/0004-637X/746/1/74
  * Roulston et al. (2021) Roulston, B. R., Green, P. J., Toonen, S., et al. 2021, arXiv:2105.00037
  * Sahai, Vlemmings, & Nyman (2017) Sahai R., Vlemmings W. H. T., Nyman L.-Å., 2017, ApJ, 841, 110. doi:10.3847/1538-4357/aa6d86
  * Sand et al. (2020) Sand C., Ohlmann S. T., Schneider F. R. N., Pakmor R., Röpke F. K., 2020, A&A, 644, A60. doi:10.1051/0004-6361/202038992
  * Schreier, Hillel & Soker (2019) Schreier R., Hillel S., Soker N., 2019, MNRAS, 490, 4748
  * Schreier & Soker (2016) Schreier R., Soker N., 2016, RAA, 16, 70. doi:10.1088/1674-4527/16/5/070
  * Shiber & Iaconi (2021) Shiber S., Iaconi R., 2021, AAS
  * Shiber et al. (2019) Shiber S., Iaconi R., De Marco O., Soker N., 2019, MNRAS, 488, 5615
  * Shiber, Kashi, & Soker (2017) Shiber S., Kashi A., Soker N., 2017, MNRAS, 465, L54. doi:10.1093/mnrasl/slw208
  * Shiber, Schreier & Soker (2016) Shiber S., Schreier R., Soker N., 2016, RAA, 16, 117
  * Shiber & Soker (2018) Shiber S., Soker N., 2018, MNRAS, 477, 2584. doi:10.1093/mnras/sty843
  * Shore et al. (1994) Shore S. N., Livio M., van den Heuvel E. P. J., Nussbaumer H., Orr A., 1994, inbi.conf
  * Soker (2004) Soker N., 2004, NewA, 9, 399. doi:10.1016/j.newast.2004.01.004
  * Soker (2015) Soker N., 2015, ApJ, 800, 114. doi:10.1088/0004-637X/800/2/114
  * Soker (2016) Soker N., 2016, NewA, 47, 16. doi:10.1016/j.newast.2016.01.005
  * Soker (2020) Soker N., 2020, Galax, 8, 26. doi:10.3390/galaxies8010026
  * Soker & Gilkis (2018) Soker N., Gilkis A., 2018, MNRAS, 475, 1198
  * Soker, Grichener & Gilkis (2019) Soker N., Grichener A., Gilkis A., 2019, MNRAS, 484, 4972
  * Staff et al. (2016) Staff J. E., De Marco O., Macdonald D., Galaviz P., Passy J.-C., Iaconi R., Low M.-M. M., 2016, MNRAS, 455, 3511. doi:10.1093/mnras/stv2548
  * Valenti, Ferraro, & Origlia (2004) Valenti E., Ferraro F. R., Origlia L., 2004, MNRAS, 351, 1204. doi:10.1111/j.1365-2966.2004.07861.x
  * van den Heuvel (1976) van den Heuvel E. P. J., 1976, IAUS, 73, 35
  * Vigna-Gómez, et al. (2020) Vigna-Gómez A., et al., 2020, arXiv, arXiv:2001.09829
  * Whitehouse et al. (2021) Whitehouse L. J., Farihi J., Howarth I. D., Mancino S., Walters N., Swan A., Wilson T. G., et al., 2021, MNRAS, 506, 4877. doi:10.1093/mnras/stab1913
  * Zou et al. (2022) Zou Y., Chamandy L., Carroll-Nellenback J., Blackman E. G., Frank A., 2022, arXiv, arXiv:2202.05715


# PauliComposer: Compute Tensor Products of Pauli Matrices Efficiently

Sebastián V. Romero  0000-0002-4675-4452<EMAIL_ADDRESS>TECNALIA, Basque Research and Technology Alliance (BRTA), 48160 Derio, Spain
Department of Physical Chemistry, University of the Basque Country UPV/EHU,
Apartado 644, 48080 Bilbao, Spain  Juan Santos-Suárez  0000-0001-9360-2411
<EMAIL_ADDRESS>Instituto Galego de Física de Altas Enerxías (IGFAE),
Universidade de Santiago de Compostela, 15705 Santiago de Compostela, Spain

(August 29, 2024)

###### Abstract

We introduce a simple algorithm that efficiently computes tensor products of
Pauli matrices. This is done by tailoring the calculations to this specific
case, which allows to avoid unnecessary calculations. The strength of this
strategy is benchmarked against state-of-the-art techniques, showing a
remarkable acceleration. As a side product, we provide an optimized method for
one key calculus in quantum simulations: the Pauli basis decomposition of
Hamiltonians.

tensor product, Kronecker product, Pauli matrices, quantum mechanics, quantum
computing.

##  I Introduction

Pauli matrices [1] are one of the most important and well-known set of
matrices within the field of quantum physics. They are particularly important
both in physics and chemistry when used to describe Hamiltonians of many-body
spin glasses [2, 3, 4, 5, 6, 7] or for quantum simulations [8, 9, 10, 11, 12,
13]. The vast majority of these systems are out of analytic control so that
they are usually simulated through exact diagonalization which requires their
Hamiltonians to be written in its matrix form. While this task may be regarded
as a trivial matter in a mathematical sense, it involves the calculation of an
exponentially growing number of operations. Furthermore, description of
quantum systems via Matrix Product States (MPS) [14], Density Matrix
Renormalization Group (DMRG) [15] and Projected Entangled Pair States (PEPS)
[16] also involve large scale Hamiltonians, as well as Lanczos method [17],
whose formulation has been efficiently encoded on quantum hardware recently
[18].

In this work, we present the _PauliComposer_ (PC) algorithm which
significantly expedites this calculation. It exploits the fact that any Pauli
word only has one element different from zero per row and column, so a number
of calculations can be avoided. Additionally, each matrix entry can be
computed without performing any multiplications. Even though the exponential
scaling of the Hilbert space cannot be avoided, PC can boost inner
calculations where several tensor products involving Pauli matrices appear. In
particular, those that appear while building Hamiltonians as weighted sums of
Pauli strings or decomposing an operator in the Pauli basis.

The PC algorithm could be implemented in computational frameworks in which
this sort of operations are crucial, such as the Python modules Qiskit [19],
PennyLane [20], OpenFermion [21] and Cirq [22]. It can also potentially be
used in many other applications, such as the Pauli basis decomposition of the
Fock space [23] and conventional computation of Ising model Hamiltonians to
solve optimization problems [24, 25, 26, 27], among others.

The rest of the article is organized as follows: in Sec. II we describe the
algorithm formulation in depth, showing a pseudocode-written routine for its
computation. In Sec. III, a set of tests is performed to show that a
remarkable speed-up can be achieved when compared to state-of-the-art
techniques. In Sec. IV, we show how this PC algorithm can be used to solve
relevant problems. Finally, the conclusions drawn from the presented results
are given in Sec. V. We provide proofs for several statements and details of
the algorithm in the appendices.

##  II Algorithm formulation

In this section we discuss the PC algorithm formulation in detail. Pauli
matrices are hermitian, involutory and unitary matrices that together with the
identity form the set $\sigma_{\\{0,1,2,3\\}}={\\{I,X,Y,Z\\}}$. Given an input
string $x=x_{n-1}\dots x_{0}\in\\{0,1,2,3\\}^{n}$, the PC algorithm constructs

$P(x)\coloneqq\sigma_{x_{n-1}}\otimes\sigma_{x_{n-2}}\otimes\dots\otimes\sigma_{x_{0}}.$
(1)

Let us denote its matrix elements as $P_{j,k}(x)$ with $j,k=0,\dots,2^{n}-1$.
It is important to remark that for each row $j$, there will be a single column
$k(j)$ such that $P_{j,k(j)}\neq 0$ (see App. A). The solution amounts to a
map from the initial Pauli string to the positions and values of the $2^{n}$
nonzero elements. This calculation will be done sequentially, hence the
complexity of the algorithm will be bounded from below by this number.

As a first step, it is worth noting that Pauli string matrices are either real
(all elements are $\pm 1$) or purely imaginary (all are $\pm i$). This depends
on $n_{Y}$, the number of $Y$ operators in $P(x)$. We can redefine
$\tilde{Y}\coloneqq iY$, so that
$\smash[t]{\tilde{\sigma}_{\\{0,1,2,3\\}}=\\{I,X,\tilde{Y},Z\\}}$ and
$\smash[t]{\tilde{P}(x)\coloneqq\tilde{\sigma}_{x_{n-1}}\otimes\dots\otimes\tilde{\sigma}_{x_{0}}}$.
As a result, every entry in $\smash[t]{\tilde{P}(x)}$ will be $\pm 1$. This
implies that there is no need to compute any multiplication: the problem
reduces to locating the nonzero entries in $\smash[t]{\tilde{P}(x)}$ and
tracking sign changes. The original $P(x)$ can be recovered as
$\smash[t]{P(x)=(-i)^{n_{Y}\text{ mod }4}\tilde{P}(x)}$.

We will now present an iterative procedure to compute $\tilde{P}$ by finding
for each row $j$ the nonzero column number $k(j)$ and its corresponding value
$\smash[t]{\tilde{P}_{j,k(j)}}$. For the first row, $j=0$, the nonzero element
$\smash[t]{\tilde{P}_{0,k(0)}}$, can be found at

$k(0)=[y(x_{n-1})\dots y(x_{0})]_{10},$ (2)

where $[\,\cdot\,]_{10}$ is the decimal representation of a bit string and
$y(x_{i})$ tracks the diagonality of $\sigma_{x_{i}}$, where $y(x_{i})$ is
equal to $0$ if $x_{i}=\\{0,3\\}$ (thus $\sigma_{x_{i}}\in\\{I,Z\\}$) and $1$
otherwise (thus $\sigma_{x_{i}}\in\\{X,Y\\}$). The value of this entry is

$\tilde{P}_{0,k(0)}=+1\implies P_{0,k(0)}=(-i)^{n_{Y}\text{ mod }4}.$ (3)

The following entries can be computed iteratively. At the end of stage $l$,
with $l=0,\cdots,n-1$, all nonzero elements in the first $2^{l+1}$ rows of
$\smash[t]{P_{j,k(j)}}$ will have been computed using the information given by
the substring $x_{l}\dots x_{0}$. At the next step, $l+1$, the following
$2^{l}$ rows are filled using the ones that had already been computed, where
the row-column relation $k(j)$ is given by

$k(j+2^{l})=k(j)+(-1)^{y(x_{l})}2^{l},\quad j=0,\dots,2^{l}-1.$ (4)

The second term of the RHS of this relation takes into account the way that
the blocks of zeros returned at stage $l$ affect the new relative location of
the nonzero blocks within the new $2^{l+1}\times 2^{l+1}$ subcomposition. Its
corresponding values are obtained from the previous ones, up to a possible
change of sign given by

$P_{j+2^{l},k(j+2^{l})}=\epsilon_{l}P_{j,k(j)},$ (5)

with $\epsilon_{l}$ equal to $1$ if $x_{l}\in\\{0,1\\}$ and $-1$ otherwise.
This $\epsilon_{l}$ is nothing but a parameter that takes into account if
$\sigma_{x_{l}}$ introduces a sign flip. In Alg. 1 a pseudocode that
summarizes the presented algorithm using (2)-(5), is shown.

For the particular case of diagonal Pauli strings (only $I$ and $Z$ matrices),
there is no need to compute the row-column relation $k(j)$, just the sign
assignment is enough. Even if this is also the case for anti-diagonal
matrices, we focus on the diagonal case due to its relevance in combinatorial
problems [24, 25, 26, 27]. See Alg. 2 for the pseudocode of this case (PDC
stands for _PauliDiagonalComposer_).

input : $x_{n-1}x_{n-2}\dots x_{0}\leftarrow$ string with
$x_{i}\in\\{0,1,2,3\\}$

1 $n\leftarrow\textnormal{{len(}}\textnormal{\emph{x}}\textnormal{{)}}$

2 $n_{Y}\leftarrow\,$number of $Y$ matrices in $x$

$j\leftarrow\textnormal{{range(}}\textnormal{\emph{$0,2^{n}-1$}}\textnormal{{)}}$

// rows

$k,m\leftarrow\,$empty $2^{n}$-array

// columns/entries

3 $k(0)\leftarrow y(x_{n-1})\dots y(x_{0})$ in base 10

4 $m(0)\leftarrow(-i)^{n_{Y}\text{ mod }4}$

5 for  _$l\in$ range(_$0,n-1$_)_ do

6 $k(2^{l}:2^{l+1}-1)\leftarrow k(0:2^{l}-1)+(-1)^{y(x_{l})}2^{l}$

7 if  _$x_{l}\in\\{0,1\\}$_ then // $\epsilon_{l}=1$

8 $m(2^{l}:2^{l+1}-1)\leftarrow m(0:2^{l}-1)$

9 else // $\epsilon_{l}=-1$

10 $m(2^{l}:2^{l+1}-1)\leftarrow-m(0:2^{l}-1)$

11

12

output : $P(x)$ as a sparse matrix stacking $(j,k,m)$

Algorithm 1 PC: compose $n$ Pauli matrices

input : $x_{n-1}x_{n-2}\dots x_{0}\leftarrow$ string with $x_{i}\in\\{0,3\\}$

1 $n\leftarrow\textnormal{{len(}}\textnormal{\emph{x}}\textnormal{{)}}$

$j,k\leftarrow\textnormal{{range(}}\textnormal{\emph{$0,2^{n}-1$}}\textnormal{{)}}$

// rows/columns

$m\leftarrow\,$empty $2^{n}$-array

// entries

2 $m(0)\leftarrow 1$

3 for  _$l\in$ range(_$0,n-1$_)_ do

4 if  _$x_{l}=0$_ then // $\epsilon_{l}=1$

5 $m(2^{l}:2^{l+1}-1)\leftarrow m(0:2^{l}-1)$

6 else // $\epsilon_{l}=-1$

7 $m(2^{l}:2^{l+1}-1)\leftarrow-m(0:2^{l}-1)$

8

9

output : $P(x)$ as a sparse matrix stacking $(j,k,m)$

Algorithm 2 PDC: compose $n$ diagonal Pauli matrices

The PC algorithm is able to circumvent the calculation of a significant amount
of operations. When generic Kronecker product routines (see App. B) are used
for the same task, the amount of multiplications needed for computing a Pauli
string is $\mathcal{O}[n2^{2n}]$ and $\mathcal{O}[n2^{n}]$ for dense and
sparse matrices, respectively. In contrast, the PC algorithm, considering the
worst-case scenarios, needs

  * •

$\\{I,Z\\}^{\otimes n}$: $\mathcal{O}[2^{n}]$ changes of sign.

  * •

Otherwise: $\mathcal{O}[2^{n}]$ sums and $\mathcal{O}[2^{n}]$ changes of sign.

In all cases our algorithm can significantly outperform those that are not
specifically designed for Pauli matrices.

On top of that, this method is also advantageous for computing weighted Pauli
strings. Following (3), $W\coloneqq\omega P$, with arbitrary $\omega$, can be
computed by defining $W_{0,k(0)}=\omega(-i)^{n_{Y}\text{ mod }4}$ which avoids
having to do any extra multiplication. This change is reflected in Alg. 1 by
changing line 1 to $m(0)\leftarrow\omega(-i)^{n_{Y}\text{ mod }4}$ and line 2
to $m(0)\leftarrow\omega$ in Alg. 2. This is specially important as it can be
used to compute Hamiltonians written as a weighted sum of Pauli strings, where
$H=\sum_{x}\omega_{x}P(x)$.

##  III Benchmarking

In this section we analyse the improvement that the PC strategy introduces
against other known algorithms labelled as Naive (regular Kronecker product),
Algorithm 993 (Alg993) [28], Mixed and Tree [29, 30]. Further details can be
found in App. B. We benchmark these algorithms using MATLAB [31] as it is
proficient at operating with matrices (it incorporates optimized routines of
the well-known BLAS and LAPACK libraries [32, 33]). The PC avoids matrix
operations and thus it would not be ideal to implement it using MATLAB.
Instead, we use Python [34] since many quantum computing libraries are written
in this language [19, 22, 21, 20]. See Tab. 1 for a full description of the
computational resources used.

Table 1: Computer specifications. Processor | Intel Core i7-11850H ($16\times$2.50\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$)
---|---
RAM |  $32.0\text{\,}\mathrm{G}\mathrm{B}$ (DDR4)
OS | Ubuntu 22.04.1 LTS ($\times$64)
MATLAB [31] | 9.12.0.1884302 (R2022a)
Python [34] | 3.9.12
NumPy [35] | 1.23.2 | SciPy [36] | 1.9.0
Qiskit [19] | 0.38.0 | PennyLane [20] | 0.23.1
Figure 1: (Color online) Execution times for computing general (solid line)
and diagonal (dashed) $n$-Pauli strings using different methods.

Concerning memory needs, with this algorithm only $2^{n}$ nonzero elements out
of $2^{2n}$ are stored. This is exactly the same as using sparse matrices,
thus, no major improvement is to be expected. As for the computational time,
we compare how different algorithms behave as the length $n$ of the Pauli
string increases. In Fig. 1 execution times for general and diagonal Pauli
strings are shown. For the PC methods, we use the PC routine (Alg. 1) for the
general case and the PDC routine (Alg. 2) for the diagonal one. In accordance
to our theoretical analysis, the PC algorithm proves to be the best performing
routine.

On a more technical note, when using the PC routine, matrices with complex
values ($n_{Y}$ odd) take twice as much time as real valued ones ($n_{Y}$
even). Consequently, we compute their execution times separately and then
average them. Moreover, it is convenient to choose when to use PC or PDC as
the latter can be up to 10 times faster.

##  IV Real use cases of the PauliComposer algorithm

The PC algorithm can be used to perform useful calculations in physics. In
this section, the Pauli basis decomposition of a Hamiltonian and the
construction of a Hamiltonian as a sum of weighted Pauli strings are discussed
in detail. Another worth mentioning scenario is the digital implementation of
the complex exponential of a Pauli string, i.e. $e^{-i\theta
P(x)}=\cos(\theta)I-i\sin(\theta)P(x)$.

_Pauli basis decomposition of a Hamiltonian_.—The decomposition of a
Hamiltonian written as a $2^{n}\times 2^{n}$ matrix into the Pauli basis is a
common problem in quantum computing. Given a general Hamiltonian $H$, this
decomposition can be written as $H=\sum_{x}\omega_{x}P(x)$ with
$x=x_{n-1}\dots x_{0}$ and $P(x)$ as in (1). The coefficients $\omega_{x}$ are
obtained from the orthogonal projection as

$\omega_{x}=\frac{1}{2^{n}}\mathrm{tr}\\!\left[P(x)H\right]=\frac{1}{2^{n}}\sum_{j=0}^{2^{n}-1}P_{j,k(j)}(x)H_{k(j),j}.$
(6)

Following the discussion in Sec. II, the double sum collapses to a single one
in (6) since there is only one nonzero element per row and column. Each of
these weights can be computed independently, which allows for a parallel
implementation. Additionally, in some special cases, it can be known in
advance if some $\omega_{x}$ will vanish:

  * •

If $H$ is symmetric, strings with an odd number of $Y$ matrices can be avoided
($2^{n-1}(2^{n}+1)$ terms).

  * •

If $H$ is diagonal, only strings composed by $I$ and $Z$ will contribute
($2^{n}$ terms).

Table 2: Execution times (in seconds) for decomposing an arbitrary $2^{n}\times 2^{n}$ matrix. Here, PC and PDC calculations were made computing weights sequentially and in parallel. $n$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---
Non-hermitian matrix $H_{\text{NH}}$
PC (sequential) | $0.0005$ | $0.0021$ | $0.012$ | $0.078$ | $0.55$ | $4.06$ | $31.2$ | $254$ | $2008$
PC (parallel) | $0.094$ | $0.093$ | $0.11$ | $0.15$ | $0.38$ | $2.10$ | $13.5$ | $94.3$ | $719$
Qiskit | $0.0015$ | $0.0050$ | $0.020$ | $0.14$ | $1.16$ | $8.78$ | $92.38$ | $1398$ | $26938$
Hermitian matrix $H_{\text{H}}$
PC (sequential) | $0.0004$ | $0.0021$ | $0.012$ | $0.078$ | $0.56$ | $4.24$ | $32.86$ | $261$ | $2007$
PC (parallel) | $0.068$ | $0.070$ | $0.079$ | $0.12$ | $0.33$ | $1.99$ | $13.02$ | $96.5$ | $647$
Qiskit | $0.0010$ | $0.0035$ | $0.018$ | $0.10$ | $1.47$ | $12.02$ | $108$ | $1295$ | $26848$
PennyLane | $0.0013$ | $0.0060$ | $0.030$ | $0.15$ | $2.23$ | $10.66$ | $97.6$ | $2019$ | $35014$
Symmetric matrix $H_{\text{S}}$
PC (sequential) | $0.0003$ | $0.0010$ | $0.0058$ | $0.036$ | $0.24$ | $1.78$ | $14.05$ | $108$ | $794$
PC (parallel) | $0.059$ | $0.059$ | $0.061$ | $0.078$ | $0.13$ | $0.48$ | $2.75$ | $20.1$ | $140$
Qiskit | $0.0010$ | $0.0036$ | $0.018$ | $0.10$ | $1.45$ | $11.07$ | $105$ | $1320$ | $26399$
PennyLane | $0.0011$ | $0.0054$ | $0.027$ | $0.13$ | $1.36$ | $9.22$ | $91.52$ | $1477$ | $31583$
Diagonal matrix $H_{\text{D}}$
PDC (sequential) | $0.0001$ | $0.0002$ | $0.0006$ | $0.0018$ | $0.0068$ | $0.025$ | $0.094$ | $0.37$ | $1.49$
PDC (parallel) | $0.055$ | $0.057$ | $0.059$ | $0.060$ | $0.060$ | $0.064$ | $0.078$ | $0.12$ | $0.35$
Qiskit | $0.0010$ | $0.0035$ | $0.018$ | $0.10$ | $1.46$ | $11.0$ | $103$ | $1270$ | $25977$
PennyLane | $0.0010$ | $0.0047$ | $0.023$ | $0.11$ | $1.20$ | $8.29$ | $86.2$ | $1370$ | $30941$
Figure 2: (Color online) Execution times for decomposing $2^{n}\times 2^{n}$
(a) non-hermitian $H_{\text{NH}}$, (b) hermitian $H_{\text{H}}$, (c) symmetric
$H_{\text{S}}$ and (d) diagonal $H_{\text{D}}$ matrices with different
methods. For PC and PDC, solid (dotted) line depicts sequential (parallelized)
decomposition. See Tab. 2. As expected, notice that the larger $n$, the higher
impact of parallelization.

The operations made by _PauliDecomposer_ (PD) are

  * •

If $H$ is diagonal ($\mathcal{O}[2^{n}]$ strings): $\mathcal{O}[2^{2n}]$
operations.

  * •

Otherwise ($\mathcal{O}[2^{2n}]$ strings): $\mathcal{O}[2^{3n}]$ operations.

This PD algorithm checks if the input matrix satisfies one of the
aforementioned cases and computes the coefficients using the PC routine and
(6), discarding all vanishing Pauli strings. This workflow considerably
enhances our results, especially for diagonal matrices.

In Tab. 2 and Fig. 2, we tested the most extended methods for decomposing
matrices into weighted sums of Pauli strings against PD, using Python [34] to
compare their performance. In particular, we used the SparsePauliOp class from
Qiskit [19] and the decompose_hamiltonian function from PennyLane [20] (only
works with hermitian Hamiltonians). To the best of authors’ knowledge, both
routines are based on Naive approach without inspecting the input matrix
nature before proceeding.

Four types of random $2^{n}\times 2^{n}$ matrices were generated, namely non-
hermitian $H_{\text{NH}}$, hermitian $H_{\text{H}}$, symmetric $H_{\text{S}}$
and diagonal $H_{\text{D}}$ matrices. The PD vastly outperforms Qiskit and
PennyLane routines, specially for the symmetric and diagonal cases.

_Building of a Hamiltonian as a sum of weighted Pauli strings_.—Many
Hamiltonians are written in terms of weighted Pauli strings. Our method can
compute weighted Pauli strings directly without extra computations. In Fig. 3
we show a performance comparison of the presented methods for computing
Hamiltonians written as sums of weighted Pauli strings. The Hamiltonian used
is similar to the one proposed in [27],

$H=\sum_{i=0}^{n-1}\alpha_{i}\sigma^{i}_{3}+\sum_{i<j}^{n-1}\beta_{ij}\sigma^{i}_{3}\sigma^{j}_{3},$
(7)

being the corresponding weights $\vec{\alpha}$ and $\vec{\beta}$ arbitrary and
$\sigma^{i}_{3}$ as defined in (11). This Hamiltonian is computed using Alg.
3, which uses the PDC routine (see Alg. 2) with two inputs: the string
$x\in\\{0,3\\}^{n}$ to compute and the weights to consider. In the PDC case,
we use two strategies: compute each weighted term of (7) directly and compute
each Pauli string and then multiply it by its corresponding weight (solid and
dashed lines in Fig. 3, respectively). This is done by changing lines 3 to
$H\leftarrow H+\alpha_{i}\texttt{PDC(}str_{1}\texttt{)}$ and 3 to $H\leftarrow
H+\beta_{ij}\texttt{PDC(}str_{2}\texttt{)}$ in Alg. 3 for the second one.
There is no significant difference between both methods.

input : $\vec{\alpha},\vec{\beta}\leftarrow\,$lists of weights

1
$n\leftarrow\textnormal{{len(}}\textnormal{\emph{$\vec{\alpha}$}}\textnormal{{)}}$

2 $H\leftarrow\,2^{n}\times 2^{n}$ sparse matrix of zeros

3 for  _$i\in$ range(_$0,n-1$_)_ do

$str_{1}\leftarrow\,$string of $n$ zeros

// $n$ identities

$str_{1}(i)\leftarrow 3$

// $Z$ in the $i$-th position

4 $H\leftarrow
H+\textnormal{{PDC(}}\textnormal{\emph{$str_{1},\alpha_{i}$}}\textnormal{{)}}$

5 for  _$j\in$ range(_$i+1,n-1$_)_ do

6
$str_{2}\leftarrow\textnormal{{copy(}}\textnormal{\emph{$str_{1}$}}\textnormal{{)}}$

$str_{2}(j)\leftarrow 3$

// $Z$ in the $j$-th position

7 $H\leftarrow
H+\textnormal{{PDC(}}\textnormal{\emph{$str_{2},\beta_{ij}$}}\textnormal{{)}}$

8

9

output : Hamiltonian $H$ as a sparse matrix

Algorithm 3 Ising model Hamiltonian computation Figure 3: (Color online)
Execution times for computing (7) using Alg. 3 (solid line) and computing
previously the Pauli string and multiply it by its corresponding weight
(dashed).

##  V Conclusions

The fast and reliable computation of tensor products of Pauli matrices is
crucial in the field of quantum mechanics and, in particular, of quantum
computing. In this article we propose a novel algorithm with proven
theoretical and experimental enhancements over similar methods of this key yet
computationally tedious task. This is achieved by taking advantage of the
properties of Pauli matrices and the tensor product definition, which implies
that one can avoid trivial operations such as multiplying constants by one and
waste time computing elements with value zero that could be known in advance.

Concerning memory resources, it is convenient to store the obtained results as
sparse matrices since only $2^{n}$ out of $2^{2n}$ entries will not be zero
for a Pauli string of length $n$, i.e. the density of the resultant matrix
will be $2^{-n}$ (see App. A).

Our benchmark tests suggest that the PauliComposer algorithm and its variants
can achieve a remarkable acceleration when compared to the most well-known
methods for the same purpose both for single Pauli strings and real use cases.
In particular, the most considerable outperformance can be seen in Tab. 2 for
the symmetric and diagonal matrix decomposition over the Pauli basis.

Finally, its simple implementation (Alg. 1-2) can potentially allow to
integrate the PC routines into quantum simulation packages to enhance inner
calculations.

###### Acknowledgements.

We thank Javier Mas Solé, Yue Ban and Mikel García de Andoin for the helpful
discussions. This research is funded by the project “BRTA QUANTUM: Hacia una
especialización armonizada en tecnologías cuánticas en BRTA” (expedient no.
KK-2022/00041). The work of JSS has received support from Xunta de Galicia
(Centro singular de investigación de Galicia accreditation 2019-2022) by
European Union ERDF, from the Spanish Research State Agency (grant
PID2020-114157GB-100) and from MICIN with funding from the European Union
NextGenerationEU (PRTR-C17.I1) and the Galician Regional Government with own
funding through the “Planes Complementarios de I+D+I con las Comunidades
Autónomas” in Quantum Communication.

Data and code availability statement. The data used in the current study is
available upon reasonable request from the corresponding authors. The code
used can be found at https://github.com/sebastianvromero/PauliComposer.

##  Appendix A Some proofs regarding Pauli strings

In this section we prove two key properties of Pauli strings on which our
algorithm is based.

######  Proposition 1.

A Pauli string $P(x)$ of length $n$ given by (1) has only $2^{n}$ nonzero
entries.

###### Proof.

With the help of Fig. 4, we can compute the number of zeros in the resulting
matrix as

$\displaystyle n_{0}(n)$
$\displaystyle={\color[rgb]{1,0,0}2\left(2^{n-1}\times
2^{n-1}\right)}+{\color[rgb]{0,0.5,0}4\left(2^{n-2}\times 2^{n-2}\right)}$ (8)
$\displaystyle+{\color[rgb]{0,0,1}8\left(2^{n-3}\times
2^{n-3}\right)}+\dots+2^{n}(1\times 1)$
$\displaystyle=\sum_{k=n}^{2n-1}2^{k}=2^{n}\left(2^{n}-1\right).$

In other words, $P(x)$ will have only $2^{n}$ nonzero terms. We can prove (8)
by induction easily: since $n_{0}(n=1)$ is true, if we assume that $n_{0}(n)$
holds we can see that

$\displaystyle n_{0}(n+1)$ $\displaystyle=2\cdot 2^{n}(2^{n}-1)+2\cdot
2^{2n}=2^{n+1}\left(2^{n+1}-1\right)$ (9)

also holds true. ∎

$\bigotimes_{i=0}^{n-1}\sigma_{x_{n-i-1}}=\begin{bmatrix}\begin{array}[]{cc}{\color[rgb]{0,0.5,0}0}&\begin{array}[]{cc}\cdots&{\color[rgb]{0,0,1}0}\\\
{\color[rgb]{0,0,1}0}&\cdots\end{array}\\\
\begin{array}[]{cc}\cdots&{\color[rgb]{0,0,1}0}\\\
{\color[rgb]{0,0,1}0}&\cdots\end{array}&{\color[rgb]{0,0.5,0}0}\\\
\end{array}&{\color[rgb]{1,0,0}0}\\\
{\color[rgb]{1,0,0}0}&\begin{array}[]{cc}{\color[rgb]{0,0.5,0}0}&\begin{array}[]{cc}\cdots&{\color[rgb]{0,0,1}0}\\\
{\color[rgb]{0,0,1}0}&\cdots\end{array}\\\
\begin{array}[]{cc}\cdots&{\color[rgb]{0,0,1}0}\\\
{\color[rgb]{0,0,1}0}&\cdots\end{array}&{\color[rgb]{0,0.5,0}0}\\\
\end{array}\\\ \end{bmatrix}$ Figure 4: (Color online) Scheme for computing
the number of zeros of an arbitrary composition of $n$ Pauli matrices.

######  Corollary 1.1.

A Pauli string $P(x)$ of length $n$ given by (1) has only one nonzero entry
per row and column.

###### Proof.

Since the tensor product of unitary matrices is also unitary, then $|\\!\det
P(x)|=1$. From Th. 1, only $2^{n}$ entries of the resulting $2^{n}\times
2^{n}$ matrix are nonzero. So the logical conclusion to be drawn is that the
unique way to locate them without having a row and a column full of zeros,
thus returning a zero determinant, is that each row and column must have only
one nonzero entry. ∎

##  Appendix B Standard methods for computing tensor products

In this appendix we briefly review the well established algorithms that were
used in the benchmark [28, 29, 30]. First, one can consider what we call the
Naive algorithm, which consists on performing the calculations directly. It is
clearly highly inefficient as it scales in the number of operations as
$\mathcal{O}[n2^{n}]$ for sparse Pauli matrices. Second, the Mixed algorithm
uses the mixed-product property

$\bigotimes_{i=0}^{n-1}\sigma_{x_{n-i-1}}=\prod_{i=0}^{n-1}\sigma^{i}_{x_{n-i-1}},$
(10)

with

$\sigma^{i}_{x_{i}}\coloneqq\begin{cases}I^{\otimes
n-1}\otimes\sigma_{x_{0}}&\text{if }i=0\\\ I^{\otimes
n-i-1}\otimes\sigma_{x_{i}}\otimes I^{\otimes i}&\text{if }0<i<n-1\\\
\sigma_{x_{n-1}}\otimes I^{\otimes n-1}&\text{if }i=n-1\end{cases},$ (11)

to simplify the calculation into a simple product of block diagonal matrices.
Based on this procedure, Alg993 is presented in [28]. It can be shown that
this method performs over $\mathcal{O}[n2^{n}]$ operations. Besides that, as
Fig. 1 suggests, the fact that it requires to transpose and reshape several
matrices has a non-negligible effect that fatally increases its computation
time. Finally, the Tree routine starts storing pairs of tensor products as

$\begin{aligned}
\displaystyle\left\\{\sigma_{x_{n-2i-1}}\otimes\sigma_{x_{n-2i-2}}\right\\}^{n/2-1}_{i=0}&&&\text{if
}n\text{ even}\\\
\displaystyle\left\\{\sigma_{x_{n-1}}\right\\}\cup\left\\{\sigma_{x_{n-2i-1}}\otimes\sigma_{x_{n-2i-2}}\right\\}^{\lfloor
n/2\rfloor}_{i=0}&&&\text{if }n\text{ odd}\end{aligned},$ (12)

and proceeds with the resultant matrices following the same logic, which
allows to compute (1) by iteratively grouping its terms by pairs. For better
results, this method can be parallelized.

## References

  * Pauli [1927] W. Pauli, Zeitschrift für Physik 43, 601 (1927).
  * Heisenberg [1928] W. Heisenberg, Zeitschrift für Physik 49, 619 (1928).
  * Bethe [1931] H. Bethe, Zeitschrift für Physik 71, 205 (1931).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett. 35, 1792 (1975).
  * Panchenko [2012] D. Panchenko, Journal of Statistical Physics 149, 362 (2012).
  * Hubbard and Flowers [1963] J. Hubbard and B. H. Flowers, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Altland and Simons [2006] A. Altland and B. Simons, Second quantization, in  _Condensed Matter Field Theory_ (Cambridge University Press, 2006) pp. 39–93.
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Zeitschrift für Physik 47, 631 (1928).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210 (2002).
  * Seeley  _et al._ [2012] J. T. Seeley, M. J. Richard, and P. J. Love, The Journal of Chemical Physics 137, 224109 (2012).
  * Tranter  _et al._ [2015] A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert, F. Wilhelm, and P. J. Love, International Journal of Quantum Chemistry 115, 1431 (2015).
  * Tranter  _et al._ [2018] A. Tranter, P. J. Love, F. Mintert, and P. V. Coveney, Journal of Chemical Theory and Computation 14, 5617 (2018).
  * Steudtner and Wehner [2018] M. Steudtner and S. Wehner, New Journal of Physics 20, 063010 (2018).
  * Östlund and Rommer [1995] S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, arxiv:cond-mat/0407066 (2004).
  * Lanczos [1950] C. Lanczos, J. Res. Natl. Bur. Stand. B 45, 255 (1950).
  * Kirby  _et al._ [2023] W. Kirby, M. Motta, and A. Mezzacapo, Quantum 7, 1018 (2023).
  * Qiskit [2021] Qiskit, Qiskit: An Open-source Framework for Quantum Computing (2021).
  * PennyLane [2018] PennyLane, PennyLane: Automatic Differentiation of Hybrid Quantum-Classical Computations (2018).
  * OpenFermion [2017] OpenFermion, OpenFermion: The Electronic Structure Package for Quantum Computers (2017).
  * Cirq [2022] Cirq, Cirq (2022).
  * Liu  _et al._ [2022] R. Liu, S. V. Romero, I. Oregi, E. Osaba, E. Villar-Rodriguez, and Y. Ban, Entropy 24, 1529 (2022).
  * Lucas [2014] A. Lucas, Frontiers in Physics 2, 5 (2014).
  * Osaba  _et al._ [2022] E. Osaba, E. Villar-Rodriguez, and I. Oregi, IEEE Access 10, 55805 (2022).
  * V. Romero  _et al._ [2023] S. V. Romero, E. Osaba, E. Villar-Rodriguez, I. Oregi, and Y. Ban, Scientific Reports 13, 11777 (2023).
  * G. de Andoin  _et al._ [2022] M. G. de Andoin, E. Osaba, I. Oregi, E. Villar-Rodriguez, and M. Sanz, in  _Proceedings of the Genetic and Evolutionary Computation Conference Companion_, GECCO ’22 (2022) pp. 2214–2222.
  * Fackler [2019] P. L. Fackler, ACM Trans. Math. Softw. 45, 1 (2019).
  * Horn and Johnson [1991] R. A. Horn and C. R. Johnson, Matrix equations and the kronecker product, in  _Topics in Matrix Analysis_ (Cambridge University Press, 1991) p. 239–297.
  * Burrus [2009] C. S. Burrus, Implementing Kronecker Products Efficiently, in  _Automatic Generation of Prime Length FFT Programs_ (OpenStax CNX, 2009) pp. 41–49.
  * MAT [2022] _MATLAB version 9.12.0.1884302 (R2022a)_, The Mathworks, Inc., Natick, Massachusetts (2022).
  * Lawson  _et al._ [1979] C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh, ACM Trans. Math. Softw. 5, 308 (1979).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK users’ guide_ , 3rd ed., Software, environments, tools (Society for Industrial and Applied Mathematics, 1999).
  * Python [2022] Python,  _Python: A Dynamic, Open Source Programming Language_, Python Software Foundation (2022).
  * Charles R. Harris and K. Jarrod Millman et al. [2020] Charles R. Harris and K. Jarrod Millman et al., Nature 585, 357 (2020).
  * SciPy [2020] SciPy, Nature Methods 17, 261 (2020).