# Are Interactions with Neutron Star Merger Winds Shaping the Jets?

L. Nativi1, G. P. Lamb2, S. Rosswog1, C. Lundman1, G. Kowal3,
1Department of Astronomy and Oskar Klein Centre, Stockholm University,
AlbaNova 10691 Stockholm, Sweden
2School of Physics and Astronomy, University of Leicester, University Road,
Leicester LE1 7RH, UK
3Escola de Artes, Ciências e Humanidades, Universidade de São Paulo, Av.
Arlindo Béttio, 1000 – Vila Guaraciaba,
CEP: 03828-000, São Paulo – SP, Brazil E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Jets can become collimated as they propagate through dense environments and
understanding such interactions is crucial for linking physical models of the
environments to observations. In this work, we use 3D special-relativistic
simulations to study how jets propagate through the environment created around
a neutron star merger remnant by neutrino-driven winds. We simulate four jets
with two different initial structures, top-hat and Gaussian, and two
luminosities. After jet breakout, we study the angular jet structures and the
resulting afterglow light curves. We find that the initial angular structures
are efficiently washed out during the propagation, despite the small wind mass
of only $\sim 10^{-3}$ M⊙. The final structure depends on the jet luminosity
as less energetic jets are more strongly collimated, and entrainment of
baryons leads to a moderate outflow Lorentz factor ($\approx 40$). Although
our jets are not specifically intended to model the outflows of the GW170817
event, we show that they can be used to produce light curves consistent with
the afterglow observed in the aftermath of GW170817. Using this procedure we
show how the inferred physical parameters e.g., inclination angle, ambient
particle number density, can vary substantially between independent fits of
the same dataset and appear to be sensitive to smaller details of the angular
jet shape, indicating that observationally inferred parameters may depend
sensitively on the employed jet models.

###### keywords:

gamma-ray bursts – method: numerical – hydrodynamics – jets and outflows –
neutron star mergers – relativistic processes

††pubyear: 2021††pagerange: Are Interactions with Neutron Star Merger Winds
Shaping the Jets?–Are Interactions with Neutron Star Merger Winds Shaping the
Jets?

##  1 Introduction

Binary neutron star (BNS) mergers have long been suspected to produce the
central engines of short gamma ray bursts (sGRB) (Eichler et al., 1989). The
link was firmly established in August 2017, after the combined detection of
gravitational waves and a short GRB from the same BNS merger (Abbott et al.,
2017a, b, c; Goldstein et al., 2017; Savchenko et al., 2017). The actual
origin of the gamma-ray signal is still debated, coming from either an off-
axis jet or from the shock breakout of a relativistic cocoon inflated by the
jet itself (Lundman & Beloborodov, 2021; Gottlieb et al., 2018a; Gottlieb et
al., 2018b; Kasliwal et al., 2017; Mooley et al., 2018a; Nakar & Piran, 2017).
Nevertheless, the multi-band observations of a rising afterglow (Margutti et
al., 2017, 2018; D’Avanzo et al., 2018; Lyman et al., 2018; Lamb et al., 2019;
Alexander et al., 2017; Hallinan et al., 2017; Troja et al., 2017, 2018, 2019)
together with the detection of superluminal motion (Mooley et al., 2018b;
Ghirlanda et al., 2019; Hotokezaka et al., 2019) have settled the presence of
a jet that successfully broke out from the surrounding ejecta, observed off-
axis with a viewing angle $\theta_{\rm obs}\approx 19^{\circ}$ (Mooley et al.,
2018a; Lazzati et al., 2018; Murguia-Berthier et al., 2017; Lamb et al., 2018;
Margutti & Chornock, 2020). The information obtainable from afterglow
observations is strongly dependent on the angular structure of the emerging
jet. This structure might be mainly determined by the launching process
(Kathirgamaraju et al., 2019) or it may arise as a consequence of the
interaction with the surrounding environment during propagation. Therefore, an
understanding of the processes that shape the jet could in principle provide
insights into both jet formation, and the post-merger environment. Recent
relativistic (magneto-) hydrodynamics simulations (Lazzati et al., 2021;
Urrutia et al., 2021; Pavan et al., 2021; Geng et al., 2019; Gottlieb et al.,
2020, 2021b; Gottlieb & Nakar, 2021; Gottlieb et al., 2021a; Nathanail et al.,
2021; Murguia-Berthier et al., 2014; Murguia-Berthier et al., 2017, 2021b;
Beniamini et al., 2020b; Hamidani & Ioka, 2021) have illustrated the
importance of the ambient medium in shaping the jet, thereby highlightning the
importance of understanding the remnant structure and ejecta properties. The
joint events GW170817 and GRB170817A were followed by an additional
electromagnetic transient spanning the spectral bands from UV to optical and
IR on time scales from days to weeks (Abbott et al., 2017b; Arcavi et al.,
2017; Cowperthwaite et al., 2017; Evans et al., 2017; Drout et al., 2017; Pian
et al., 2017; Smartt et al., 2017; Soares-Santos et al., 2017; Tanvir et al.,
2017; Utsumi et al., 2017). The observed properties were consistent with the
expectations for a thermal transient powered by the radioactive decay of
freshly synthesized r-process elements (a so-called "macronova" or "kilonova"
e.g., Li & Paczyński, 1998; Kulkarni, 2005; Rosswog, 2005; Metzger et al.,
2010; Roberts et al., 2011; Kasen et al., 2013; Kasen et al., 2015, 2017; Yu
et al., 2013; Metzger, 2017; Tanaka et al., 2017; Perego et al., 2017; Rosswog
et al., 2018). Understanding the properties of this signal requires in-depth
investigation of all the processes that can unbind material during and after a
BNS merger, together with the available amount of free neutrons provided by
each ejection channel. By now, several mass-ejection channels have been
identified and they differ in terms of launch time, mass, electron fraction
$Y_{e}$ and velocity. During the merger $\sim 10^{-3}-10^{-2}$ M⊙ of material
are ejected dynamically (Rosswog et al., 1998, 1999; Rosswog & Davies, 2002;
Oechslin et al., 2007; Bauswein et al., 2013; Radice et al., 2018a). On longer
time scales ($\sim 1$ s) a few $10^{-2}$ M⊙ of material can be unbound from
the torus surrounding the remnant by the action of nuclear heating (Fernandez
& Metzger, 2013; Fernandez et al., 2015; Just et al., 2015; Metzger et al.,
2008), magnetic (Siegel & Ciolfi, 2015; Siegel & Metzger, 2017, 2018; Ciolfi
et al., 2017; Ciolfi & Kalinani, 2020; Murguia-Berthier et al., 2021a) and
viscous effects (Shibata et al., 2017; Shibata & Hotokezaka, 2019; Radice et
al., 2018b; Fujibayashi et al., 2018, 2020). Weak interactions play a key role
since they can change the initially extremely low electron fraction and they
are therefore of paramount importance for nucleosynthesis and the
electromagnetic appearance of a BNS merger (Ruffert et al., 1997; Rosswog &
Liebendörfer, 2003; Dessart et al., 2009; Perego et al., 2014, 2017; Martin et
al., 2015, 2018; Miller et al., 2019; Murguia-Berthier et al., 2021a).

In Nativi et al. 2021 we investigated how the propagation of a jet through a
neutrino-driven wind (Perego et al., 2014) impacts the resulting radioactive
transient. We found in particular that an outer layer of low-$Y_{e}$ material
with only very small mass can have a significant impact on the resulting
electromagnetic signal, once more underlining how crucial the realistic
modelling of the jet environment is for interpreting observations (see also
Lazzati et al., 2021; Pavan et al., 2021; Ito et al., 2021). For example, a
relativistic jet can punch away such low-$Y_{e}$, high-opacity material
causing the macronova to become brighter earlier for observers close to the
polar axis.
Here, we investigate the same neutrino-driven wind environment as in Nativi et
al. 2021, based on the simulations of Perego et al. (2014), but from a
different perspective. We want to know to which extent the initial jet
structure is preserved when a jet propagates through a realistic post-merger
environment. To this end we perform high-resolution 3D special relativistic
hydrodynamic simulations of jets with different initial angular structures and
luminosities. In Section 2 we briefly discuss the numerical methods adopted
for this study, and the simulation setup. We present our findings in Section 3
and, as a conclusive exercise, we attempt a comparison between the light
curves produced by our resulting jet structures and the only afterglow with a
certain connection to neutron star mergers i.e., GRB 170817A. We summarize and
list our main findings in Section 4.

##  2 Methods

We use the same special-relativistic hydrodynamic high-order shock-capturing
adaptive mesh refinement Godunov-type code AMUN111AMUN code is open source and
freely available from https://gitlab.com/gkowal/amun-code., as in Nativi et
al. 2021. In the present simulations the Riemann states at each cell interface
are reconstructed using the 5th-order Monotonicity Preserving (MP5) scheme of
Suresh & Huynh (1997), and the Harten-Lax-van Leer (HLL) Riemann solver
(Harten et al., 1983) is used. The time integration is performed using the 3rd
order Strong Stability Preserving Runge-Kutta (SSPRK) algorithm from Gottlieb
et al. (2011) with a value of 0.5 for the Courant-Friedrichs-Lewy number
(Courant et al., 1967). The numerical integration of the special relativistic
hydrodynamic equations requires nonlinear conversion of the conservative to
primitive variables. In this work we used the $1D_{W}$ scheme with a Newton-
Raphson iterative solver (see, e.g., Noble et al., 2006) with the tolerance of
$10^{-8}$ as the primitive variable solver.

###  2.1 Simulations setup

We explore here how a jet is impacted by the interaction with the ambient
medium produced by a neutrino-driven wind. We use the wind profiles of Perego
et al. (2014) in all simulations. We embed the wind in a steeply decreasing
power-law atmosphere that settles to a constant density plateau with
$\rho_{\rm atm}=10^{-12}$ g cm3. The background density has been chosen low
enough to ensure that the final bow-shock produced by the jet carries a
negligible amount of mass-energy. This requirement has been later verified at
the end of each simulation.

The computational grid is shaped as a rectangular Cartesian box above the
equatorial plane, with extensions $-4\times 10^{4}\leq(x\rm{,}y)\leq 4\times
10^{4}$ km in the polar direction and $0\leq z\leq 1.6\times 10^{5}$ km in the
axial direction. We used the properties of our mesh to set up a grid
characterized by a superposition of two patches. One patch is kept fixed from
the beginning and is used to cover the early jet propagation inside the wind,
keeping the maximum available resolution ($2.4\times 2.4\times 3.2$ km). This
patch has a radial extension of $R_{\rm p1}=100$ km and an axial extension of
$z_{\rm p1}=2000$ km. The second patch has an initial resolution of four times
that of the first patch and a radial extension222 We use the adaptive mesh
refinement algorithm included in AMUN to ensure that the difference in
resolution between neighbouring computational blocks never exceeds a factor of
2. $R_{\rm p2}=2000$ km. In the axial direction the second patch is used to
dynamically follow the jet head defining $z_{\rm p2}=(t-t_{0})c$ for all
$z_{\rm p2}>z_{\rm p1}$, with $c$ the speed of light and $t_{0}$ the actual
time at which the jet is launched.

####  2.1.1 Jets

As in Nativi et al. 2021, we inject the jet as a boundary condition from the
lower boundary of the domain (Gottlieb et al., 2021b, 2020, 2018a; Harrison et
al., 2018; Geng et al., 2019; Mizuta & Ioka, 2013; Mizuta & Aloy, 2009).
Around the origin we move the lower boundary above the equatorial plane to a
height of $z_{\rm 0}=40$ km, where injection occurs. In this way we remove the
presence of the remnant, whose evolution is not of interest in the current
work, but we keep within the domain the material around the rest of the
equatorial plane. To inject the jet into the computational domain we specify
five functions: density $\rho_{\rm j}(\theta)$ and pressure $p_{\rm
j}(\theta)$ in the local rest-frame and the three components of the velocity
$\mathbf{v}_{\rm j}(\theta)=(v_{\rm j}^{\rm x},v_{\rm j}^{\rm y},v_{\rm
j}^{\rm z})$ in the laboratory frame, $\theta$ is the angle from the jet axis.
The geometry of the jet-injection region is characterized by a spherical
radius from the source $r_{\rm 0}$ and an opening angle $\theta_{\rm j}$, with
corresponding solid angle is $\Omega_{\rm j}=2\pi(1-\cos\theta_{\rm j})$. The
intrinsic properties of the plasma flowing into the domain can be obtained
from its luminosity $L_{\rm j}$, its velocity $v_{\rm 0}^{\rm r}$ with $r$
expressing the radial direction (or equivalently the Lorentz factor
$\Gamma_{\rm 0}^{\rm r}$) and a description of the contribution of the
internal energy $h$. In this paper we present the results of four simulations,
characterized by two luminosities ($L_{\rm j}=$1050 and 1051 erg s-1) and two
different initial structures: Gaussian (gs51 and gs50) and top-hat (th51 and
th50). To set the geometry, we chose a height for the jet injection $z_{\rm
0}$ and an initial opening angle $\theta_{\rm j}=15^{\circ}$, consistent with
the results from GRMHD simulations (Kathirgamaraju et al., 2019). In all the
simulations the engine is active for $\Delta t_{\rm inj}=100$ ms with constant
luminosity, and decays exponentially after being shut-off on a time scale of
10 ms. The injected plasma at the top of the injection region is defined by an
initial radial Lorentz factor $\Gamma_{\rm 0}=5$ and a specific enthalpy in
the jet core $h_{\rm c}=30$, resulting in an asymptotic Lorentz factor
$\Gamma_{\infty}=150$. All the simulations adopt a polytropic equation of
state with adiabatic exponent 4/3, appropriate for relativistic gases.

For the initial jet structure we start from the assumption that the luminosity
per unit solid angle

$\frac{dL}{d\Omega}=r^{2}v_{\rm 0}\Gamma_{\rm 0}^{2}\rho c^{2}h,$ (1)

is the same for an observer along the axis for any structure. We call
$r=z_{0}/\cos\theta$ the injection radius i.e., the distance from the origin.
We proceed with a few simplifying assumptions: the local speed is radial,
$v_{0}\neq v_{0}(\theta)$, $\rho_{\rm j}\neq\rho_{\rm j}(\theta)$ (for fixed
$r$) while the enthalpy

$h=1+\frac{4p_{\rm j}}{\rho_{\rm j}},$ (2)

is allowed to have an angular dependence. For a top-hat jet, $h=h_{\rm c}$ is
constant and for a Gaussian jet, $h(\theta=0)=h_{\rm c}$ i.e., the Gaussian
jet has the same $h$ along its axis as the top-hat jet. In this approach we
are assuming that the whole inner structure can be ascribed to the thermal
energy of the jet. The total jet luminosity is then

$L_{\rm j}=\int\frac{dL}{d\Omega}d\Omega=2\pi\int\frac{dL}{d\Omega}d\mu,$ (3)

where $\mu\equiv\cos\theta$.

For a top-hat jet the luminosity per unit solid angle is constant for
$\theta<\theta_{j}$, and the total luminosity is

$L=2\pi\int\frac{dL}{d\Omega}d\mu=2\pi(1-\mu_{j})\frac{dL}{d\Omega},$ (4)

where $\mu_{j}=\cos\theta_{j}$. We can then solve for $\rho_{\rm j}$ at a
given $r$ using Equation 1,

$\rho_{\rm j}=\frac{L_{\rm
j}}{2\pi(1-\mu_{j})r^{2}v_{0}\Gamma_{0}^{2}h_{c}c^{2}},$ (5)

and find the pressure, $p_{\rm j}$, from Equation 2,

$p_{\rm j}=\frac{1}{4}(h_{\rm c}-1)\rho_{\rm j}c^{2}.$ (6)

The primitive variables $(v,\rho,p)$ are thus fully specified for the top-hat
jet.

For the Gaussian case we start again with our assumptions $v_{0}\neq
v_{0}(\theta)$ and $\rho_{\rm j}\neq\rho_{\rm j}(\theta)$. Therefore, using
the same equation for $\rho$ (Equation 5), we only need to find the new
equation for $p=p(\theta)$ that corresponds to a Gaussian jet with the same
pressure along the axis as the top-hat jet, and also carries the same total
luminosity.

We want $h$ to have a Gaussian profile, therefore we assume

$h(\theta)-1=(h_{\rm
c}-1)\exp\left(-\frac{\theta^{2}}{\theta_{j}^{2}}\right).$ (7)

So, to find the primitive variables for a Gaussian jet, $v=v_{0}$, $\rho$ is
computed from Equation 5 (same as for the top-hat jet) and $p$ is found by
taking $h$ from Equation 7 and subbing into Equation 6. The top-hat and
Gaussian jet shapes are shown in Figure 4.

##  3 Results

Overall, we find the dynamical evolution of the system consistent with the
expectations from other results (Mizuta & Aloy, 2009; Bromberg et al., 2011;
Mizuta & Ioka, 2013; Nagakura et al., 2014; Murguia-Berthier et al., 2014;
Murguia-Berthier et al., 2021b; Harrison et al., 2018; Gottlieb et al., 2018a,
2020, 2021b; Gottlieb & Nakar, 2021; Duffell et al., 2018; Nativi et al.,
2021; Urrutia et al., 2021; Ito et al., 2021; Hamidani et al., 2020; Hamidani
& Ioka, 2021). In all models, the jets are quickly engulfed by a cocoon and
experience several recollimation shocks before they manage to break out from
the ejecta. Figures 1, 2 and 3 show the evolution of all models for rest-mass
density, Lorentz factor and $(h-1)$ respectively. After breaking out, all the
previously collimated jets start expanding sideways widening a funnel with an
opening angle close to the initial $\theta_{j}$. Such angular extension is
also visible in the final angular profiles and is consistent with the
expectations (e.g. Mizuta & Ioka, 2013; Nagakura et al., 2014; Harrison et
al., 2018), see Fig. 4. The post-breakout evolution proceeds with a
progressive sideways expansion of the jet accompanied by a gradual decreasing
of $(h-1)$ and an increasing Lorenz factor $\Gamma$. We run all simulations
until the jet reaches a configuration of ballistic expansion, with $(h-1)<1$
everywhere on the grid. After this point, essentially all the internal energy
in the plasma has been converted to kinetic energy, and the system undergoes
free expansion. For all models we choose to extract the angular profiles at
$\approx 0.4$ s.

In order to use the final jet structure to produce afterglow light curves, we
need to estimate the properties of the forward afterglow shock. Specifically,
we need the energy per solid angle d$E(\theta)$/d$\Omega$, and the Lorentz
factor of the shock as a function of angle $\overline{\Gamma}(\theta)$. At
late times, the whole jet has been compressed into a thin layer behind the
forward shock. Therefore, we estimate these properties by radially averaging
the relevant quantities, while keeping their angular dependences. We use mass-
averages in order to conserve both mass and energy in the averaging process.

However, the averaging process should not consider all plasma inside the grid.
For instance, slow, heavy plasma that flows behind the jet will not have time
to catch up and contribute to the shock on relevant time scales. We therefore
need a procedure to estimate which parts of the jet are likely to contribute
to the afterglow shock at relevant times. Here, we simply adopt a threshold
value of $\Gamma h=2$; fluid with a lower value is assumed to not contribute
much to the afterglow333Afterglow lightcurves calculated using a threshold
$\Gamma h=5$, for the Lorentz factor profile, show no significant deviation
when compared to those found using our adopted limit for fits to GW170817
data, see §3.1. Low-$\Gamma$ and low-energy components, that constitute the
wider cocoon, will contribute to the afterglow for a mildly inclined system at
the afterglow onset time, and for observers at much higher system inclinations
will become the dominant emission component.. This value ensures that the
mass-averaged Lorentz factor is not affected much by the slow plasma that
follows the jet, while retaining plasma at fairly large angles to the jet
axis. We perform a cylindrical average of the 3D grid returned by the
hydrodynamic simulation: this reduces the problem to 2D in the $(r,z)$ plane
and with uniform spacing $\Delta r$, $\Delta z$. The physical quantities we
need are:

$D=\Gamma\rho\;\mbox{,}\quad h=1+\frac{4P}{\rho c^{2}}\;\mbox{,}\quad
E=D\Gamma hc^{2}-p-Dc^{2}$ (8)

at every grid point. We assume the polar angle $\theta$ being distributed with
uniform spacing $\Delta\theta$ in the interval $[0,\pi/2]$, then iterate over
the grid. As motivated above, we adopt a threshold value of $\Gamma h=2$ to
allow the selection of relativistic material from the jet and cocoon, which
may contribute to the final afterglow emission. For each selected cell we then
compute the corresponding polar angle for the angular bin
$[\theta_{l};\theta_{l+1}]$. Where for each angular bin we define two
variables:

$\displaystyle\Delta E_{l}=\Delta E_{l}+E_{ij}\mathcal{V}_{ij}$ (9)
$\displaystyle\Delta M_{l}=\Delta M_{l}+D_{ij}\mathcal{V}_{ij}$ (10)
$\displaystyle\Delta\Gamma_{l}=\Delta\Gamma_{l}+\Gamma_{ij}D_{ij}\mathcal{V}_{ij}$
(11)

where $\mathcal{V}_{ij}=2\pi r_{i}\Delta r\Delta z$ is the volume of cell
$\\{i,j\\}$. Finally, we can get the two required quantities as:

$\displaystyle\frac{\rm{d}E(\theta)}{\rm{d}\Omega}\approx\frac{\Delta
E_{l}}{\Delta\Omega_{l}}=\frac{\Delta E_{l}}{2\pi\lvert\Delta\mu\rvert}\
\rm{,with}\quad\mu=\cos(\theta),$ (12)
$\displaystyle\overline{\Gamma}(\theta)\approx\frac{\Delta\Gamma_{l}}{\Delta
M_{l}}.$ (13)

The resulting profiles are shown in Fig 4. Clearly, there are structural
differences, but in all cases the initial jet structure is largely washed out
in the interaction with the ejecta despite their relatively small amount of
mass ($\sim 10^{-3}$ M⊙).

Figure 1: Vertical slices ($y=0$) of the rest-mass density shortly after the
breakout (left) and after the jet shut-off (right) for all models. In all
cases the jet punches a hole into the surrounding wind roughly proportional to
its luminosity, consistently with what shown in Nativi et al. 2021. In the
lateral direction the wind keeps expanding homologously. Particularly in the
two high luminosity models (upper figures) some low density material is
preceded by some at higher density. The material sitting on the top is
injected before the breakout and has experienced more mixing. The matter
injected after the breakout is only marginally affected by the interaction and
maintains a lower density. These effects appear being more important in the
Gaussian model, and are also visible in Fig. 2 and Fig 3. Figure 2: Same as in
Fig 1 but for the Lorentz factor. All jets experience a first strong
collimation shock converting the specific kinetic energy into thermal and then
several recollimation events (see Fig 3). After the breakout the jet drives a
strong relativistic shock in the surrounding environment. At later times the
jet material keeps being adiabatically accelerated towards its asymptotic
value $\Gamma h$. The material injected before the breakout is affected by
baryon pollution and is characterised by an averagely lower value for
$\Gamma$, as described in Fig 1. Figure 3: Same as in Fig. 1 and Fig. 2 but
showing the specific enthalpy $h-1$. These plots together with those of the
Lorentz factor describe the evolution of the jet. It is clear how the most
polluted material is characterised by both a lower $\Gamma$ and a lower $h-1$,
meaning a sensitive reduction of $\Gamma h$ due to mass entrainment for the
material injected before the breakout.

The Gaussian jets have a resultant profile with a broader low energy wing and
they are less energetic on-axis, while the top-hat jets result in a profile
that is slightly more energetic and faster at small angles suggesting that for
such initial jet structures it is easier to preserve the energetic spine. In
all models we find that, due to baryon entrainment, the final average Lorentz
factor is systematically lower than their asymptotic value
($\Gamma_{\infty}=\Gamma h\approx 150$). In the late time maps, shown on the
right of Fig. 3, the material just below the head is characterised by a much
lower specific enthalpy than the material injected later. All the matter that
is flowing into the domain before the jet breakout is heavily affected by the
interaction with its surroundings. This interaction triggers hydrodynamic
instabilities and turbulence at the jet-ejecta interface, favouring mixing and
mass entrainment. The material injected after the jet breakout is, however,
only marginally affected by the instabilities at the interface, and this
material proceeds along the opened up funnel retaining most of its energy.

###  3.1 Application to the GRB 170817A afterglow

Figure 4: The resultant energy (solid red line) and Lorentz factor (solid blue
line) profiles with polar angle for each model at time t $\approx 0.4$s. The
top four panels include the injected energy (dashed red) and Lorentz factor
(dotted blue) profiles for each model. The bottom panels show the ratio of the
resultant energies with angle between an initially Gaussian profile (middle
panels) and an initially ‘top-hat’ profile (top panels) jet. For each profile,
$\theta_{j}=0.25\sim 15^{\circ}$, marked in the bottom panel by a vertical
dotted line, the horizontal dashed line indicates gs/th$=1$. The resultant
energy and Lorentz factor profiles are calculated by summing cells with
$\Gamma h\geq 2$. The jet opening angle as inferred from the jet-break timing
for on-axis afterglow light curves is shown for each model as a vertical dash-
dotted grey line.

The resulting energy and Lorentz factor jet profiles will be effectively
frozen after break out and the jet becomes conical, expanding as a function of
distance from the central engine. This resultant jet structure can be used to
inform the choice of fiducial energy and Lorentz factor profiles for
structured jet afterglow models, such as those used to model the late time
afterglow to GRB 170817A (e.g. Lyman et al., 2018; Troja et al., 2018, 2019;
Lamb & Kobayashi, 2018; Lamb et al., 2019; Resmi et al., 2018; Ryan et al.,
2020; Salafia et al., 2019; Fernández et al., 2021).

We use the energy and Lorentz factor profiles shown in Fig. 4 with the
afterglow models developed in Lamb & Kobayashi (2017); Lamb et al. (2018);
Lamb et al. (2021) to predict and model the appearance of GRB afterglow
lightcurves from structured jets at any inclination. Although it is not our
intention to infer the afterglow properties of GRB 170817A, we use this
event’s afterglow light curve to demonstrate that our resultant jet profiles
can nevertheless produce afterglows consistent with the expectation from
structured jet models. We use two cases for each model, without lateral
spreading, where the jet is assumed to have no sideways expansion i.e.,
perpendicular to the radial direction, and with lateral spreading, where we
assume sideways expansion at the local sound-speed of the jet element (see
Lamb et al., 2018; Lamb et al., 2021, for details). We use the X-ray, optical,
and radio frequency data-sets as presented in Troja et al. (2021);
Balasubramanian et al. (2021); Makhathini et al. (2020) and references
therein444Note that in all cases, the very late time X-ray and 3 GHz radio
data, as presented here, can be accommodated by the GRB afterglow model with
either no, or moderate, lateral spreading – where sound-speed lateral
spreading is assumed to represent an upper limit on the sideways expansion..
As the energy and initial Lorentz factors are set by our resultant profiles,
we fix the electron index, $p=2.15$, following standard closure relations for
the afterglow spectrum, and adopt a redshift, $z=0.009783$ (Levan et al.,
2017) with Planck cosmology giving a distance, $D_{L}\sim 43.7$ Mpc. We fit
for the inclination, $\iota$, the ambient particle number density, $n$, and
the microphysical parameters, $\varepsilon_{B}$ and $\varepsilon_{e}$. The
best fit parameters and lightcurves are found via a Markov Chain Monte Carlo
(MCMC) using emcee (Foreman-Mackey et al., 2013), with flat log priors for
each parameter except the inclination, which has a flat prior in $\cos\iota$.
The prior ranges and the central plus 16th and 84th percentile values are
shown in Table 1. A sample of lightcurves drawn randomly from the posterior
distribution are shown in Fig. 5 with the observed data points, where the
models without lateral spreading have solid lines, and those with maximal
lateral spreading, have dash-dotted lines.

Table 1: Free parameters for the MCMC model lightcurve fits to GW170817/GRB 170817A afterglow data. All priors are flat in the range shown. Each model profile is fit with no lateral spreading, and with sound-speed lateral spreading. Note that the model jet structures are found via simulations which are inconsistent with the conditions in GW170817/GRB 170817A, and as such, the fit afterglow parameter values should not be used to directly model or infer the properties of this event. Model | Lat. spread (y/n) |  $\iota$ (rad) | $\log\varepsilon_{B}$ | $\log\varepsilon_{e}$ |  $\log n$ $\log$(cm-3)
---|---|---|---|---|---
Prior | – | $[\cos(0.75),~{}\cos(0.25)]$ | $[-6,~{}-0.3]$ | $[-6,~{}-0.3]$ | $[-6,~{}2]$
th50 | n | $0.29^{+0.04}_{-0.03}$ | $-2.31^{+0.62}_{-0.74}$ | $-1.54^{+0.52}_{-0.42}$ | $-3.12^{+0.49}_{-0.34}$
th50 | y | $0.25^{+0.01}_{-0.00}$ | $-0.59^{+0.05}_{-0.07}$ | $-2.20^{+0.04}_{-0.03}$ | $-4.03^{+0.06}_{-0.03}$
gs50 | n | $0.41^{+0.05}_{-0.05}$ | $-3.24^{+0.82}_{-0.65}$ | $-0.93^{+0.49}_{-0.62}$ | $-2.17^{+0.35}_{-0.44}$
gs50 | y | $0.28^{+0.03}_{-0.02}$ | $-1.01^{+0.33}_{-1.09}$ | $-1.89^{+0.82}_{-0.16}$ | $-3.96^{+0.38}_{-0.21}$
th51 | n | $0.45^{+0.02}_{-0.02}$ | $-4.37^{+1.14}_{-1.10}$ | $-1.60^{+0.75}_{-0.79}$ | $-1.36^{+0.20}_{-0.19}$
th51 | y | $0.37^{+0.02}_{-0.01}$ | $-2.78^{+0.96}_{-1.48}$ | $-1.85^{+1.06}_{-0.86}$ | $-2.66^{+0.15}_{-0.17}$
gs51 | n | $0.51^{+0.03}_{-0.02}$ | $-4.29^{+1.09}_{-1.15}$ | $-1.65^{+0.70}_{-0.76}$ | $-1.01^{+0.18}_{-0.18}$
gs51 | y | $0.43^{+0.01}_{-0.01}$ | $-3.47^{+1.03}_{-1.23}$ | $-1.73^{+0.83}_{-0.74}$ | $-2.15^{+0.13}_{-0.14}$
Figure 5: 100 randomly selected afterglow lightcurves from the posterior
distributions of MCMC fits to GRB 170817A data for each jet profile where 50
lightcurves include the effects of sound speed lateral spreading (dash-
dotted), and 50 are shown without lateral spreading (solid). The energy and
Lorentz factor profiles are fixed to those from our simulations and the
inclination, ambient density and microphysical parameters are allowed to vary.

For each of our models we use an afterglow lightcurve without lateral
spreading, and with $\iota=0$ to measure the jet break time as seen by a
distant observer. By fitting a simple smoothly broken power law function for
the model afterglow light curve, as in Lamb et al. (2021), we find the best-
fit observed jet-break time in each case, we then use the known energy in the
jet and the fixed ambient number density to estimate the effective core angle
of each profile. The inferred $\theta_{c}$ angles are shown on Fig. 4 as
vertical grey dash-dotted lines in each panel. The core values in each case
are: $\theta_{c}=[0.042,~{}0.045,~{}0.136,~{}0.170]$ rad for th50, gs50, th51,
gs51 respectively. These jet core angles as estimated from the ‘observed’ jet
break time imply that, for a fixed merger environment density into which a jet
is launched and for a fixed initial jet profile, that the opening/core angle
of the resultant jet depends on the initial jet power. We note that the peak
energy for the resultant profiles, in each case, is at a comparable level,
$\sim$few $10^{52}$ erg, in terms of isotropic equivalent energy.

Using the inferred core angles for each profile, we find the ratio
$\iota/\theta_{c}$ for the posterior distribution of inclination angles in
each case.

As can be seen in Table 1, where lateral spreading is included, the preferred
inclination of the system is reduced and as $\theta_{c}=$constant in each
case, then spreading models have a typically smaller ratio $\iota/\theta_{c}$
than the non-spreading models. Where superluminal motion is included in
combination with afterglow modelling for GRB 170817A, Nakar & Piran (2021)
found that the ratio of the inclination, $\iota$, to the jet core angle,
$\theta_{c}$, should fall in the range $4\lesssim\iota/\theta_{c}\lesssim 6$,
to be consistent with all the observations. We did not include the radio
imaging constraints in our MCMC and, as such, none of our models fits
comfortably in this range. However, for the low- and high-energy models at
$L_{\rm j}=10^{50-51}$ erg s-1, the jet structure has little impact on the
ratio $\iota/\theta_{c}$, which depends more significantly on the initial
power of the jet, and whether lateral spreading is included or not. Where
lateral spreading is included, the $\iota/\theta_{c}$ values are marginally
smaller in all cases, consistent with the findings of Fernández et al. (2021).
It is, however, not the scope of this paper to consistently fit the afterglow
to GRB 170817A. We additionally note that, despite the differences in the
resultant profiles for both fixed injected energy and fixed initial jet
structure, the afterglow lightcurves can all reproduce the observed temporal
behaviour of the late afterglow to GW170817/GRB 170817A. However, we note that
although all of these resultant jet profiles can produce afterglow lightcurves
that fit the data, and all have nearly consistent ratios, $\iota/\theta_{c}$,
for fixed jet energies, the value of the other fit parameters varies
significantly between the models. For GRB afterglow models, there are some
well known parameter degeneracies and as we employ fixed energy/Lorentz factor
profiles for each model, where all other physics is kept the same, the
variation in preferred parameters between the different structure profiles is
likely a result of these e.g., the deceleration timescale is proportional to
the ratio $(E/n)^{1/3}$, and the inferred inclination angle can depend on the
jet core size or structure profile, for more details see Beniamini et al.
(2020a); Ryan et al. (2020); Nakar & Piran (2021) and, where intrinsically
more diverse jet structure profiles are considered, Takahashi & Ioka (2021).
For a jet with $L_{\rm j}=10^{50}$ erg s-1, where the initial profile is a
top-hat jet, the inferred inclination for the data is $\iota\sim 0.25\pm 0.01$
rad ($14.3\pm 0.6$ degrees) whereas for an initially Gaussian jet structure,
$\iota\sim 0.28\pm 0.03$ ($16.0\pm 1.8$ degrees) showing that, where we use
the resultant angular profiles of realistic simulations following jet
propagation through a neutron star merger environment, then the assumed
injected jet structure can influence the inferred parameters from the
subsequent afterglow fitting. This has implications for Hubble parameter
estimates from GW counterparts where simulations and/or afterglow model fits
are used (e.g. Hotokezaka et al., 2019; Wang & Giannios, 2021; Mastrogiovanni
et al., 2021).

##  4 Conclusions

We have investigated to which extent a jet from a BNS merger maintains its
initial structure when plowing through the matter previously ejected through a
neutrino-driven wind, and the consequences for the observed afterglow
emission. To achieve this, we have performed a set of 3D special-relativistic
hydrodynamic simulations where two different profiles for the initial jet
structure were adopted, and two different luminosities. We propagated the jets
in a previously simulated neutrino-driven wind environment (which was also
used in Nativi et al. (2021)). The final angular jet profiles were recovered
after the jets broke out of the ejecta, and the profiles were used to compute
afterglow light curves. Our study was not designed to model GW170817,
nevertheless, since this is so far the only observed multi-messenger event
from a merging neutron star, it is tempting to ask if the systems we have
simulated can produce an afterglow lightcurve that is consistent with that
from GW170817/GRB 170817A. Interestingly, we find that all of our resultant
jet structures can reproduce the observed, multi-band afterglow light curve of
this event.
Our main findings are:

  * •

Despite the relatively small amount of material in such winds ($\sim 10^{-3}$
M⊙) the emerging jet appears to be entirely shaped by the interaction. Hence
jets initially injected with a different structure possess a similar shape at
the end of the simulations.

  * •

Mixing is important for the final Lorentz factor profiles. Baryon entrainment
is a consequence of hydrodynamic instabilities and turbulence arising at the
contact interface between the jet and ejecta, and cannot be properly
represented in 2D simulations (Duffell et al., 2018; Harrison et al., 2018;
Urrutia et al., 2021). We find that mixing is important for the part of the
jet injected before the jet breakout, while the latter part of the jet is
barely affected, retaining its initial composition and energy per baryon (see
Harrison et al., 2018; Gottlieb et al., 2020, 2021b; Gottlieb & Nakar, 2021,
for more details).

  * •

We find that the Lorentz factor, $\Gamma$, for all our jet models peaks at a
value consistent with or above the theoretical minimum for the short GRB
population, $\Gamma\approx 40$ (Nakar, 2007).

  * •

We further explore the jet core angles, $\theta_{\rm c}$, that would be
inferred (via the jet break time) given an observation of the afterglow. These
angles show a dependency on the injected jet power. This is a consequence of
the different degree of collimation: low energy jets are more effectively
collimated by the interaction with the environment than high energy jets.

  * •

The final jet profiles differ in their small scale features, which are likely
sensitive to simulation details such as the method of jet launching, numerical
resolution as well as the properties of the surrounding environment. Smaller
features of the jet structure can therefore not be constrained by current jet
simulations, but they significantly change our best fit values of the
afterglow parameters. Where single simulation results are used to model the
afterglow (i.e. to get precise inclination angles for use in cosmology as in
e.g., Hotokezaka et al. 2019; Wang & Giannios 2021), the small differences
arising from the choice of the initial conditions and subsequent evolution
impact the parameter values inferred by the best fit.

## Acknowledgements

This work has been supported by the Swedish Research Council (VR) under grant
number 2016- 03657_3, by the Swedish National Space Board under grant number
Dnr. 107/16, the research environment grant “Gravitational Radiation and
Electromagnetic Astrophysical Transients (GREAT)" funded by the Swedish
Research council (VR) under Dnr 2016-06012 and by the Knut and Alice
Wallenberg Foundation under Dnr KAW 2019.0112. We gratefully acknowledge
stimulating interactions from COST Action CA16104 “Gravitational waves, black
holes and fundamental physics" (GWverse) and from COST Action CA16214 “The
multi-messenger physics and astrophysics of neutron stars" (PHAROS).
We acknowledge support from the Swedish National Space Agency.
The simulations were performed on resources provided by the Swedish National
Infrastructure for Computing (SNIC) at Beskow and Tetralith and on the
resources provided by the North-German Supercomputing Alliance (HLRN).
G.P.L. is supported by the Science and Technology Facilities Council, UK via
grant ST/S000453/1.
G.K. acknowledges support from the São Paulo Research Foundation, FAPESP
(grants 2013/10559-5 and 2019/03301-8).
The authors thank the anonymous referee for their timely and constructive
comments that have improved the paper.

## Data Availability

The data used to produce the observational findings, i.e. the angular profiles
describing the jets structures, are available from the author, L.N., upon
reasonable request.

## References

  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, ApJ, 848, L12
  * Abbott et al. (2017c) Abbott B. P., et al., 2017c, ApJ, 848, L13
  * Alexander et al. (2017) Alexander K. D., et al., 2017, ApJ, 848, L21
  * Arcavi et al. (2017) Arcavi I., et al., 2017, Nature, 551, 64
  * Balasubramanian et al. (2021) Balasubramanian A., et al., 2021, ApJL, 914, L20
  * Bauswein et al. (2013) Bauswein A., Goriely S., Janka H. T., 2013, ApJ, 773, 78
  * Beniamini et al. (2020a) Beniamini P., Granot J., Gill R., 2020a, MNRAS, 493, 3521
  * Beniamini et al. (2020b) Beniamini P., Duran R. B., Petropoulou M., Giannios D., 2020b, ApJ, 895, L33
  * Bromberg et al. (2011) Bromberg O., Nakar E., Piran T., Sari R., 2011, ApJ, 740, 100
  * Ciolfi & Kalinani (2020) Ciolfi R., Kalinani J. V., 2020, ApJ, 900, L35
  * Ciolfi et al. (2017) Ciolfi R., Kastaun W., Giacomazzo B., Endrizzi A., Siegel D. M., Perna R., 2017, Phys. Rev. D, 95, 063016
  * Courant et al. (1967) Courant R., Friedrichs K., Lewy H., 1967, IBM Journal of Research and Development, 11, 215
  * Cowperthwaite et al. (2017) Cowperthwaite P. S., et al., 2017, ApJ, 848, L17
  * D’Avanzo et al. (2018) D’Avanzo P., et al., 2018, A&A, 613, L1
  * Dessart et al. (2009) Dessart L., Ott C. D., Burrows A., Rosswog S., Livne E., 2009, ApJ, 690, 1681
  * Drout et al. (2017) Drout M. R., et al., 2017, Science, 358, 1570
  * Duffell et al. (2018) Duffell P. C., Quataert E., Kasen D., Klion H., 2018, ApJ, 866, 3
  * Eichler et al. (1989) Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126
  * Evans et al. (2017) Evans P. A., et al., 2017, Science, 358, 1565
  * Fernandez & Metzger (2013) Fernandez R., Metzger B. D., 2013, MNRAS, 435, 502
  * Fernandez et al. (2015) Fernandez R., Quataert E., Schwab J., Kasen D., Rosswog S., 2015, MNRAS, 449, 390
  * Fernández et al. (2021) Fernández J. J., Kobayashi S., Lamb G. P., 2021, arXiv e-prints, p. arXiv:2101.05138
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fujibayashi et al. (2018) Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, ApJ, 860, 64
  * Fujibayashi et al. (2020) Fujibayashi S., Shibata M., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., 2020, Phys. Rev. D, 101, 083029
  * Geng et al. (2019) Geng J.-J., Zhang B., Kölligan A., Kuiper R., Huang Y.-F., 2019, ApJ, 877, L40
  * Ghirlanda et al. (2019) Ghirlanda G., et al., 2019, Science, 363, 968
  * Goldstein et al. (2017) Goldstein A., et al., 2017, ApJ, 848, L14
  * Gottlieb & Nakar (2021) Gottlieb O., Nakar E., 2021, arXiv e-prints, p. arXiv:2106.03860
  * Gottlieb et al. (2011) Gottlieb S., Ketcheson D., Shu C.-W., 2011, Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations. WORLD SCIENTIFIC (https://www.worldscientific.com/doi/pdf/10.1142/7498), doi:10.1142/7498, https://www.worldscientific.com/doi/abs/10.1142/7498
  * Gottlieb et al. (2018a) Gottlieb O., Nakar E., Piran T., 2018a, MNRAS, 473, 576
  * Gottlieb et al. (2018b) Gottlieb O., Nakar E., Piran T., Hotokezaka K., 2018b, MNRAS, 479, 588
  * Gottlieb et al. (2020) Gottlieb O., Levinson A., Nakar E., 2020, MNRAS, 495, 570
  * Gottlieb et al. (2021a) Gottlieb O., Lalakos A., Bromberg O., Liska M., Tchekhovskoy A., 2021a, arXiv e-prints, p. arXiv:2109.14619
  * Gottlieb et al. (2021b) Gottlieb O., Nakar E., Bromberg O., 2021b, MNRAS, 500, 3511
  * Hallinan et al. (2017) Hallinan G., et al., 2017, Science, 358, 1579
  * Hamidani & Ioka (2021) Hamidani H., Ioka K., 2021, MNRAS, 500, 627
  * Hamidani et al. (2020) Hamidani H., Kiuchi K., Ioka K., 2020, MNRAS, 491, 3192
  * Harrison et al. (2018) Harrison R., Gottlieb O., Nakar E., 2018, MNRAS, 477, 2128
  * Harten et al. (1983) Harten A., Lax P. D., Leer B. v., 1983, SIAM Review, 25, 35
  * Hotokezaka et al. (2019) Hotokezaka K., Nakar E., Gottlieb O., Nissanke S., Masuda K., Hallinan G., Mooley K. P., Deller A. T., 2019, Nature Astronomy, 3, 940
  * Ito et al. (2021) Ito H., Just O., Takei Y., Nagataki S., 2021, arXiv e-prints, p. arXiv:2105.09323
  * Just et al. (2015) Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H.-T., 2015, MNRAS, 448, 541
  * Kasen et al. (2013) Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25
  * Kasen et al. (2015) Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777
  * Kasen et al. (2017) Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature, 551, 80
  * Kasliwal et al. (2017) Kasliwal M. M., et al., 2017, Science, 358, 1559
  * Kathirgamaraju et al. (2019) Kathirgamaraju A., Tchekhovskoy A., Giannios D., Barniol Duran R., 2019, MNRAS, 484, L98
  * Kulkarni (2005) Kulkarni S. R., 2005, arXiv e-prints, pp astro–ph/0510256
  * Lamb & Kobayashi (2017) Lamb G. P., Kobayashi S., 2017, MNRAS, 472, 4953
  * Lamb & Kobayashi (2018) Lamb G. P., Kobayashi S., 2018, MNRAS, 478, 733
  * Lamb et al. (2018) Lamb G. P., Mandel I., Resmi L., 2018, MNRAS, 481, 2581
  * Lamb et al. (2019) Lamb G. P., et al., 2019, ApJ, 870, L15
  * Lamb et al. (2021) Lamb G. P., Kann D. A., Fernández J. J., Mandel I., Levan A. J., Tanvir N. R., 2021, MNRAS
  * Lazzati et al. (2018) Lazzati D., Perna R., Morsony B. J., Lopez-Camara D., Cantiello M., Ciolfi R., Giacomazzo B., Workman J. C., 2018, Phys. Rev. Lett., 120, 241103
  * Lazzati et al. (2021) Lazzati D., Perna R., Ciolfi R., Giacomazzo B., Lopez-Camara D., Morsony B., 2021, arXiv e-prints, p. arXiv:2107.08053
  * Levan et al. (2017) Levan A. J., et al., 2017, ApJ, 848, L28
  * Li & Paczyński (1998) Li L.-X., Paczyński B., 1998, ApJ, 507, L59
  * Lundman & Beloborodov (2021) Lundman C., Beloborodov A. M., 2021, ApJ, 907, L13
  * Lyman et al. (2018) Lyman J. D., et al., 2018, Nature Astronomy, 2, 751
  * Makhathini et al. (2020) Makhathini S., et al., 2020, arXiv e-prints, p. arXiv:2006.02382
  * Margutti & Chornock (2020) Margutti R., Chornock R., 2020, arXiv e-prints, p. arXiv:2012.04810
  * Margutti et al. (2017) Margutti R., et al., 2017, ApJ, 848, L20
  * Margutti et al. (2018) Margutti R., et al., 2018, ApJ, 856, L18
  * Martin et al. (2015) Martin D., Perego A., Arcones A., Thielemann F. K., Korobkin O., Rosswog S., 2015, ApJ, 813, 2
  * Martin et al. (2018) Martin D., Perego A., Kastaun W., Arcones A., 2018, Classical and Quantum Gravity, 35, 034001
  * Mastrogiovanni et al. (2021) Mastrogiovanni S., Duque R., Chassande-Mottin E., Daigne F., Mochkovitch R., 2021, A&A, 652, 13
  * Metzger (2017) Metzger B. D., 2017, Living Reviews in Relativity, 20, 3
  * Metzger et al. (2008) Metzger B. D., Piro A. L., Quataert E., 2008, MNRAS, 390, 781
  * Metzger et al. (2010) Metzger B. D., et al., 2010, MNRAS, 406, 2650
  * Miller et al. (2019) Miller J. M., et al., 2019, Phys. Rev. D, 100, 023008
  * Mizuta & Aloy (2009) Mizuta A., Aloy M. A., 2009, ApJ, 699, 1261
  * Mizuta & Ioka (2013) Mizuta A., Ioka K., 2013, ApJ, 777, 162
  * Mooley et al. (2018a) Mooley K. P., et al., 2018a, Nature, 554, 207
  * Mooley et al. (2018b) Mooley K. P., et al., 2018b, Nature, 561, 355
  * Murguia-Berthier et al. (2014) Murguia-Berthier A., Montes G., Ramirez-Ruiz E., De Colle F., Lee W. H., 2014, ApJ, 788, L8
  * Murguia-Berthier et al. (2017) Murguia-Berthier A., et al., 2017, ApJ, 835, L34
  * Murguia-Berthier et al. (2021a) Murguia-Berthier A., et al., 2021a, arXiv e-prints, p. arXiv:2106.05356
  * Murguia-Berthier et al. (2021b) Murguia-Berthier A., Ramirez-Ruiz E., De Colle F., Janiuk A., Rosswog S., Lee W. H., 2021b, ApJ, 908, 152
  * Nagakura et al. (2014) Nagakura H., Hotokezaka K., Sekiguchi Y., Shibata M., Ioka K., 2014, ApJ, 784, L28
  * Nakar (2007) Nakar E., 2007, Phys. Rep., 442, 166
  * Nakar & Piran (2017) Nakar E., Piran T., 2017, ApJ, 834, 28
  * Nakar & Piran (2021) Nakar E., Piran T., 2021, ApJ, 909, 114
  * Nathanail et al. (2021) Nathanail A., Gill R., Porth O., Fromm C. M., Rezzolla L., 2021, MNRAS, 502, 1843
  * Nativi et al. (2021) Nativi L., Bulla M., Rosswog S., Lundman C., Kowal G., Gizzi D., Lamb G. P., Perego A., 2021, MNRAS, 500, 1772
  * Noble et al. (2006) Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641, 626
  * Oechslin et al. (2007) Oechslin R., Janka H., Marek A., 2007, A & A, 467, 395
  * Pavan et al. (2021) Pavan A., Ciolfi R., Vijay Kalinani J., Mignone A., 2021, arXiv e-prints, p. arXiv:2104.12410
  * Perego et al. (2014) Perego A., Rosswog S., Cabezón R. M., Korobkin O., Käppeli R., Arcones A., Liebendörfer M., 2014, MNRAS, 443, 3134
  * Perego et al. (2017) Perego A., Radice D., Bernuzzi S., 2017, ApJ, 850, L37
  * Pian et al. (2017) Pian E., et al., 2017, Nature, 551, 67
  * Radice et al. (2018a) Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018a, ApJ, 869, 130
  * Radice et al. (2018b) Radice D., Perego A., Hotokezaka K., Bernuzzi S., Fromm S. A., Roberts L. F., 2018b, ApJL, 869, L35
  * Resmi et al. (2018) Resmi L., et al., 2018, ApJ, 867, 57
  * Roberts et al. (2011) Roberts L. F., Kasen D., Lee W. H., Ramirez-Ruiz E., 2011, ApJ, 736, L21
  * Rosswog (2005) Rosswog S., 2005, ApJ, 634, 1202
  * Rosswog & Davies (2002) Rosswog S., Davies M. B., 2002, MNRAS, 334, 481
  * Rosswog & Liebendörfer (2003) Rosswog S., Liebendörfer M., 2003, MNRAS, 342, 673
  * Rosswog et al. (1998) Rosswog S., Thielemann F. K., Davies M. B., Benz W., Piran T., 1998, in Hillebrandt W., Muller E., eds, Nuclear Astrophysics. p. 103 (arXiv:astro-ph/9804332)
  * Rosswog et al. (1999) Rosswog S., Liebendörfer M., Thielemann F.-K., Davies M., Benz W., Piran T., 1999, A & A, 341, 499
  * Rosswog et al. (2018) Rosswog S., Sollerman J., Feindt U., Goobar A., Korobkin O., Wollaeger R., Fremling C., Kasliwal M. M., 2018, A&A, 615, A132
  * Ruffert et al. (1997) Ruffert M., Janka H., Takahashi K., Schaefer G., 1997, A & A, 319, 122
  * Ryan et al. (2020) Ryan G., van Eerten H., Piro L., Troja E., 2020, ApJ, 896, 166
  * Salafia et al. (2019) Salafia O. S., Ghirlanda G., Ascenzi S., Ghisellini G., 2019, A&A, 628, A18
  * Savchenko et al. (2017) Savchenko V., et al., 2017, ApJL, 848, L15
  * Shibata & Hotokezaka (2019) Shibata M., Hotokezaka K., 2019, Annual Review of Nuclear and Particle Science, 69, annurev
  * Shibata et al. (2017) Shibata M., Kiuchi K., Sekiguchi Y.-i., 2017, Phys. Rev. D, 95, 083005
  * Siegel & Ciolfi (2015) Siegel D. M., Ciolfi R., 2015, arXiv e-prints, p. arXiv:1505.01423
  * Siegel & Metzger (2017) Siegel D. M., Metzger B. D., 2017, Physical Review Letters, 119, 231102
  * Siegel & Metzger (2018) Siegel D. M., Metzger B. D., 2018, ApJ, 858, 52
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Soares-Santos et al. (2017) Soares-Santos M., et al., 2017, ApJ, 848, L16
  * Suresh & Huynh (1997) Suresh A., Huynh H. T., 1997, Journal of Computational Physics, 136, 83
  * Takahashi & Ioka (2021) Takahashi K., Ioka K., 2021, MNRAS, 501, 5746
  * Tanaka et al. (2017) Tanaka M., et al., 2017, PASJ, 69, 102
  * Tanvir et al. (2017) Tanvir N. R., et al., 2017, ApJ, 848, L27
  * Troja et al. (2017) Troja E., et al., 2017, Nature, 551, 71
  * Troja et al. (2018) Troja E., et al., 2018, MNRAS, 478, L18
  * Troja et al. (2019) Troja E., et al., 2019, MNRAS, 489, 1919
  * Troja et al. (2021) Troja E., et al., 2021, arXiv e-prints, p. arXiv:2104.13378
  * Urrutia et al. (2021) Urrutia G., De Colle F., Murguia-Berthier A., Ramirez-Ruiz E., 2021, MNRAS,
  * Utsumi et al. (2017) Utsumi Y., et al., 2017, PASJ, 69, 101
  * Wang & Giannios (2021) Wang H., Giannios D., 2021, ApJ, 908, 200
  * Yu et al. (2013) Yu Y.-W., Zhang B., Gao H., 2013, ApJ, 776, L40


# High-temperature self-energy corrections to x-ray absorption spectra

Tun S. Tan<EMAIL_ADDRESS>Quantum Theory Project, Department of Physics,
University of Florida, Gainesville, Florida 32611, USA.  J. J. Kas  Department
of Physics, University of Washington, Seattle, Washington 98195-1560, USA.  S.
B. Trickey  Quantum Theory Project, Department of Physics, University of
Florida, Gainesville, Florida 32611, USA.  J. J. Rehr  Department of Physics,
University of Washington, Seattle, Washington 98195-1560, USA.

(August 12, 2022)

###### Abstract

Effects of finite-temperature quasiparticle self-energy corrections to x-ray
absorption spectra are investigated within the finite-temperature
quasiparticle local density GW approximation up to temperatures $T$ of order
the Fermi temperature. To facilitate the calculations, we parametrize the
quasiparticle self-energy using low-order polynomial fits. We show that
temperature-driven decrease in the electron lifetime substantially broadens
the spectra in the near-edge region with increasing $T$. However, the
quasiparticle shift is most strongly modified near the onset of plasmon
excitations.

##  I Introduction

Calculations of x-ray absorption spectra (XAS) at finite temperature (FT) have
been carried out routinely in recent years. These studies range from
relatively low temperatures up to a few hundred K (LT), to the warm dense
matter (WDM) regime at high-temperatures (HT), where $T$ is of order the Fermi
temperature $T_{F}$ [1, 2, 3, 4, 5, 6, 7, 8, 9]. A standard FT approach is to
apply Fermi’s golden rule with initial and final states calculated using
conventional density functional theory (DFT) and Fermi occupation factors.
Since DFT is a ground state theory, FT quasiparticle corrections to DFT are
essential for HT excited state calculations [10, 11]. However, even at extreme
temperatures, e.g., many thousands of K, self-consistent field calculations of
the core-hole state have sometimes been done with ground state exchange-
correlation functionals $\varepsilon_{xc}[\rho]$ [2, 12, 13, 8]. This ground-
state approximation can be unreliable in that its validity depends strongly on
the system state and its properties. Some properties are only weakly sensitive
to the temperature dependence of exchange and correlation, at least for low
temperatures well below the WDM regime. Others, such as the electrical
conductivity and x-ray absorption spectra (XAS), are strongly temperature
dependent [14, 15]. Nevertheless, the use of temperature-dependent free-energy
exchange-correlation functionals $f_{xc}[\rho,T]$ [16, 17, 18, 19, 20, 21, 15]
alone ignores effects like inelastic losses. For example, the electron
inelastic scattering effect, which results in an energy-dependent broadening,
is often included via a post-processing step by convolution of the absorption
cross-section using an empirical model (e.g. Seah-Dench formalism [22, 23]),
or the imaginary part of the self-energy [24]. However, the finite temperature
dependence of the energy-dependent broadening is typically neglected.

Another common approach for XAS calculations has been the use of the real-
space multiple scattering (RSMS) method, which is also referred to as the
real-space Green’s function (RSGF) method [25]. This approach is the real-
space analog of the Korringa-Kohn-Rostoker (KKR) approach [26, 27, 28, 29,
30]. The method treats excited quasiparticle states via an energy-dependent
self-energy, and also takes into account the dynamic response of the system to
the suddenly created core-hole. The self-energy can be viewed as an energy-
dependent, non-local analog of the exchange-correlation potential in DFT [31].
The FT generalization of the self-energy can be done formally via the
Matsubara formalism. For example, Benedict  _et al._ used the approach to
investigate the effect of $T$ on the spectral function in jellium and
aluminum, e.g., on optical properties of solid-density Al. Alternatively, as
discussed by Kas et al. [33], the FT self-energy can be calculated using a
generalization of the Migdal approximation [34], analogous to the FT $GW$
approximation of Hedin [35].

Our main goal in this work is to discuss the effects of the FT GW self-energy
on XAS, an approach that heretofore has not been explored in detail [36]. In
particular, to facilitate the calculations we introduce a parametrization of
the quasiparticle FT GW self-energy within the $G_{0}W_{0}$ scheme [37]. As
illustrations we apply the approach to the XAS for two systems with $T$ up to
10 eV (i.e. $T$ of order $10^{5}$ K). Our calculations demonstrate that
thermal broadening due to the imaginary part of the self-energy are
significant above $T\approx 1$ eV. Although lattice vibrations also are
strongly temperature dependent, that behavior is dependent on the lattice
temperature $T_{l}$ which can differ from the electronic temperature $T$ in
non-equilibrium states, as discussed in a previous work [36].

The remainder of the paper is organized as follows. Section II. provides an
brief overview of the real-space Green’s function approach to XAS and its
dependence on the self-energy $\Sigma$. In Section III, we highlight the FT
corrections to XAS with a few examples and in Section IV, we present a brief
summary and conclusions. Throughout we use Hartree atomic units
$q_{e}=\hbar=m=1$, with $q_{e}=e$ the electron charge. Thus energies are in
Hartree and distances in Bohr, unless otherwise noted. For temperature we use
either K or eV, with 1 eV $\approx$ 11,604 K. Electron densities are expressed
in Wigner-Seitz radii $r_{s}=(3/4\pi\rho)^{1/3}$.

##  II Theory Summary

###  II.1 Finite-temperature X-ray Absorption

Formally the zero temperature X-ray absorption cross section is defined via
Fermi’s golden rule as

$\displaystyle\sigma(\omega)=4\pi^{2}\frac{\omega}{c}\sum_{i,f}|\langle\Psi_{i}|\hat{\xi}\cdot\mathbf{R}|\Psi_{f}\rangle|^{2}\delta_{\Gamma}(\omega+E_{i}-E_{f}),$
(1)

where $|\Psi_{i}\rangle$ and $|\Psi_{f}\rangle$ are the many-body initial and
final states, $\hat{\xi}$ is the polarization of the incident photon, and
$\mathbf{R}$ is the many-body position operator. Then within the single-
particle (quasiparticle) approximation with dipole interactions and the sudden
approximation, the zero temperature XAS becomes

$\displaystyle\sigma_{s}(\omega)=4\pi^{2}\frac{\omega}{c}\sum_{i,f}|\langle
i|{d}|f\rangle|^{2}\delta_{\Gamma}(\omega+\varepsilon_{i}-\varepsilon_{f}),$
(2)

where $\varepsilon_{i}$ and $\varepsilon_{f}$ are the energies of the
quasiparticle initial $|i\rangle$ and final $|f\rangle$ levels and many-body
shake-up factors $S_{0}^{2}\approx 1$ are ignored. The $\delta_{\Gamma}$
factor denotes a Lorentzian of width $\Gamma$ which includes both
quasiparticle and core-hole lifetime broadening. Here, the transition operator
${d=\hat{\xi}\cdot{\bf r}}=$ is the single-particle electric dipole operator.
The one-particle states $|i\rangle$ and $|f\rangle$ can be obtained from
Hartree-Fock theory or Kohn-Sham DFT. For the treatment via DFT, see e.g.,
Refs. 38, 39, 40.

For x-ray absorption, the number of final states $|f\rangle$ required to
compute the dipole matrix element has an impact on the computational
efficiency of evaluating Eq. (2). The present work uses the RSMS approach to
alleviate this bottleneck. In RSMS, we replace the summation over the final
states $|f\rangle$ with the retarded single-electron Green’s function
$G(\omega)$ in a basis of local site-angular momentum states $|Lj\rangle$
[25],

$\displaystyle G^{jj^{\prime}}_{LL^{\prime}}(\omega)=\sum_{f}\frac{\langle
Lj|f\rangle\langle
f|L^{\prime}j^{\prime}\rangle}{\omega-\varepsilon_{f}+i\eta}\;.$ (3)

In this expression, $j$ is the index of a given site $\mathbf{R}_{j}$ and
$L=(l,m)$ are the angular momentum quantum numbers. The initial states
$|i\rangle$ are calculated with the ground state Hamiltonian $H=p^{2}/2+v(r)$
while the final states $|f\rangle$ are described by the quasiparticle
Hamiltonian $H^{\prime}=p^{2}/2+v_{f}(r)+\Sigma(r,E)$, where $v(r)$ is the
self-consistent one-electron Hartree potential, $v_{f}$ is the final state one
electron Hartree potential in the presence of a screened core hole, and
$\Sigma$ is the dynamically screened quasiparticle self-energy discussed in
detail below.

The imaginary part of the quasiparticle self-energy $\Sigma$ accounts for the
mean free path of the photoelectron. Within the quasiparticle local density
approximation (QPLDA) [41], the self-energy is given by [42]

$\displaystyle\Sigma(\textbf{r},E,T=0)$ $\displaystyle=$ $\displaystyle
v^{LDA}_{xc}(\rho(\textbf{r}))+\Sigma_{GW}(\rho(\textbf{r}),E,T=0)$ (4)
$\displaystyle-\Sigma_{GW}(\rho(\textbf{r}),E_{F},T=0)\;.$

Here $\Sigma_{GW}$ is the $GW$ self-energy calculated at the $G_{0}W_{0}$
level of refinement, that is, without self-consistent iteration of $G$ or $W$
[35]. For simplicity from here onward, we drop the spatial dependence r.

For the FT generalization, we replace the $T=0$ $GW$ self-energy with the
finite-temperature $GW$ self-energy, $\Sigma_{GW}(T)$, and introduce
$T$-dependent Fermi occupation numbers,
$f(\varepsilon)=1/[\exp{\\{\beta(\varepsilon-\mu)\\}}+1]$ for the initial and
final states in Eq. (2). In addition, the ground state exchange-correlation
potential $v^{LDA}_{xc}$ is replaced by its FT generalization
$v^{LDA}_{xc}(T)$. Thus, the finite-temperature QPLDA self-energy is

$\displaystyle\Sigma(E,T)$ $\displaystyle=$ $\displaystyle
v^{LDA}_{xc}(\rho,T)+\Sigma_{GW}(\rho,E,T)$ (5)
$\displaystyle-\textnormal{Re}\ [\Sigma_{GW}(\rho,\mu_{T},T)]$

Lastly, by using $G^{jj^{\prime}}_{LL^{\prime}}(\omega)$ in Eq. (2) in place
of the sum over final states $|f\rangle$, the FT quasiparticle cross section
can be re-expressed as:[43]

$\displaystyle\sigma_{qp}(\omega)$ $\displaystyle=$
$\displaystyle-4\pi^{2}\frac{\omega}{c}{\textnormal{Im}}\sum_{iLL^{\prime}}\langle
i|\hat{d}G^{00}_{LL^{\prime}}(\omega+\varepsilon_{i})\hat{d}^{\dagger}|i\rangle$
(6) $\displaystyle\times
f(\varepsilon_{i})\big{[}1-f(\omega+\varepsilon_{i})\big{]},$

Here, we denote the absorbing atom by the index 0.

###  II.2 Finite-temperature Self-energy $\Sigma$

The finite-$T$ quasiparticle electron self-energy within the $GW$
approximation is defined formally [34, 44] by the expression

$\displaystyle\Sigma^{M}_{GW}({\mathbf{k}},i\omega_{m})$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta}\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\sum_{n=-\infty}^{\infty}G_{0}^{M}({\mathbf{k}}-{\mathbf{q}},i\omega_{m}-i\nu_{n})$
(7) $\displaystyle\times W^{M}({\mathbf{q}},i\nu_{n})\;.$

Here $G_{0}^{M}$ is the one-electron Matsubara Green’s function,
$W^{M}=\epsilon^{-1}v$ is the screened Coulomb interaction, and
$\omega_{m}=2(m+1)\pi k_{B}T$, $\nu_{n}=2n\pi k_{B}T$ are the Matsubara
frequencies, where $\epsilon$ is the dielectric function and $v$ is the bare
Coulomb potential. The screened interaction $W^{M}$ can be expressed in terms
of its spectral representation as

$\displaystyle W^{M}({\mathbf{q}},i\nu_{n})$ $\displaystyle=$ $\displaystyle
v({\mathbf{q}})+\int_{-\infty}^{\infty}d\omega^{\prime}\frac{D({\mathbf{q}},\omega^{\prime})}{i\nu_{n}-\omega^{\prime}+i\eta\,{\textnormal{sgn}}(\omega^{\prime})}$
(8)

where $v({\mathbf{q}})=4\pi/q^{2}$ is the bare Coulomb potential in Fourier
representation, and
$D({\mathbf{q}},\omega)=(1/\pi)|{\textnormal{Im}}\,W^{M}_{c}({\mathbf{q}},\omega)|\,{\textnormal{sgn}}(\omega)$
is the anti-symmetric (in frequency) bosonic excitation spectrum.
$W^{M}_{c}=W^{M}-v$ is the correlation part of the screened interaction.

Our choice of the electron gas dielectric function reflects a balance between
the level of physics included and computational feasibility. Thus for
simplicity, we use the random phase approximation (RPA), which is analogous to
the FT generalization of the Lindhard function [45],

$\displaystyle\epsilon({\mathbf{q}},\omega,T)$ $\displaystyle=$ $\displaystyle
1+2v({\mathbf{q}})\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\frac{f(\varepsilon_{{\mathbf{k}}-{\mathbf{q}}})-f(\varepsilon)}{\omega-\varepsilon_{{\mathbf{k}}-{\mathbf{q}}}+\varepsilon_{\mathbf{k}}+i\eta},$
(9)

where $f(\varepsilon)=1/[\exp{\\{\beta(\varepsilon-\mu)\\}}+1]$ is the Fermi-
Dirac occupation factor, and $\mu=\mu(T)$ is the chemical potential. The real
part of $\epsilon({\mathbf{q}},\omega,T)$ is obtained from the imaginary part
via a Kramers-Kronig transform. From an analytic continuation to the
real-$\omega$ axis, the FT GW retarded self-energy $\Sigma_{GW}$ is given by
the Migdal approximation [34]

$\displaystyle\Sigma_{GW}({\mathbf{k}},\omega,T)=\Sigma_{X}({\mathbf{k}},\omega,T)+\int_{0}^{\infty}d\omega^{\prime}\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}D({\mathbf{q}},\omega^{\prime})$
(10) $\displaystyle\times$
$\displaystyle\bigg{[}\frac{f(\varepsilon_{{\mathbf{k}}-{\mathbf{q}}})+N(\omega^{\prime})}{\omega+\omega^{\prime}-\varepsilon_{{\mathbf{k}}-{\mathbf{q}}}+i\eta}+\frac{1-f(\varepsilon_{{\mathbf{k}}-{\mathbf{q}}})+N(\omega^{\prime})}{\omega-\omega^{\prime}-\varepsilon_{{\mathbf{k}}-{\mathbf{q}}}+i\eta}\bigg{]},$

where
$\Sigma_{X}({\mathbf{k}},\omega,T)=\int[d^{3}{\mathbf{q}}/(2\pi)^{3}]f(\varepsilon_{{\mathbf{k}}-{\mathbf{q}}})v_{\mathbf{q}}$
is the exchange part of self-energy and
$N(\omega)=1/[\exp\\{\beta\omega\\}-1]$ is the Bose factor. The poles of the
Green’s function $G^{M}$ contribute to the Fermi occupations whereas the poles
of the screened interaction $W^{M}$ contribute to the Bose factor.

Calculations of the imaginary part of $\Sigma_{GW}({\mathbf{k}},\omega,T)$
involves a singe integral over the magnitude of ${\mathbf{q}}$ but to obtained
the real part we need to perform a Kramers-Kronig transform resulting in a
double integral. In typical RSGF XAS calculations [46], tens of thousands of
self-energy evaluations are required. Thus that calculation of quasiparticle
self-energy $\Sigma_{GW}(k,k^{2}/2,T)$ immediately becomes a computational
bottleneck. To circumvent that difficulty, we model $\Sigma_{GW}(k,k^{2}/2,T)$
via low-order polynomial fits to numerical calculation of $\Sigma_{GW}$ based
on Eq. (10), on a grid up to $T=2T_{F}$, where $T_{F}=E_{F}/k_{B}$. The form
of the fitting functions is described in detail in the Appendices. A
comparison between the QPLDA $\Sigma_{GW}$ and the fits is shown in Fig. 1 for
the homogeneous electron gas. For simplicity, we approximate the self-energy
for levels $k<k_{F}$ with $v_{xc}(T)$ independent of $k$.

Figure 1: The real part of the $\Sigma_{GW}$ (dashed) and our parametrization
(solid) for various densities: $r_{s}$ = 1.0, 2.0, 3.0, 4.0 at $T/T_{F}$ =
0.01 (blue) and 1.0 (red). For reference the LDA-$v_{xc}$(T) [20] values are
denoted as dots at $k_{F}$.

##  III Results and Discussion

In this section, we present our results showing effects of the use of the FT
$\Sigma_{GW}(T>0)$ instead of the zero $T$ values $\Sigma_{GW}(T=0)$. Note
that $\Sigma_{GW}(T=0)$ depends implicitly on the electronic temperature $T$
through the self-consistent electron density at that temperature. The FT SCF
calculations are carried out using an extension of the original RSMS code,
which is now implemented in FEFF10 [46, 47]. For the FT exchange correlation-
potential, we use the KSDT tabulation $v^{LDA}_{xc}(T)$ [20, 48] for both
calculations. In the FEFF10 calculations, we compute the atomic part with the
ground state exchange potential, whereas we use $\Sigma_{GW}(T)$ for the fine
structure.

As noted in Ref. [2, 8], corrections to the core energy levels are needed for
$T\gtrsim 1$ eV. We compute the core level shifts using the all-electron full-
potential linearized plane wave code FLEUR [49, 50, 51]. Within these
calculations, we ignore the explicit temperature dependence in the LDA
exchange correlation functional [52]. As a side note, FLEUR uses non-
overlapping muffin-tin potentials whereas FEFF uses overlapping muffin-tin
potentials.

In the RSMS formalism, the decomposition of the Green’s function $G$ into a
central atom $G^{c}$ contribution and a multiple-scattering $G^{sc}$
contribution allows us to describe the XAS in terms of the atomic background
$\sigma_{0}$ and the oscillatory fine structure $\chi$, i.e.,
$\sigma=\sigma_{0}(1+\chi)$. For many XANES calculations, it is found that the
atomic background matches the experimental results better when calculated
without the self-energy corrections due to the overestimation of the exchange
within the muffin-tin (MT) potential approximation [53]. Fig. 2 shows the
effect of using the ground-state potential versus the use of $\Sigma_{GW}(T)$
for the atomic background. The pre-edge is dominated by the atomic background
and thus is sensitive to the choice of exchange potential. The pre-edge
amplitude is reduced by $\approx 28\%$ due to the temperature correction of
$\Sigma_{GW}(T)$. More pump-probe experimental XAS measurements for $T\approx
10$ eV are required to validate the FT self-energy effect for the atomic
background. Nonetheless, the FT self-energy corrections are important for the
description of HT fine-structure.

Figure 2: $\textnormal{L}_{3}$-edge XAS for Cu (lattice constant $a=3.61$
Å[54]) at electronic temperature $T=0$ eV and 2 eV. The solid curves denote
the absorption with ground state atomic background while the dots represent
the absorption with $\Sigma_{GW}(T>0)$ self-energy atomic background. The
experimental measurement at ambient condition is shown as black cross.[55]

As a first example, we consider the FT K-edge x-ray absorption near-edge
spectrum (XANES) for aluminum (fcc Al, lattice constant $a=4.05$ Å[54]). Al is
a prototypical nearly-free electron system in the sense that the electronic
density of states (DOS) in the conduction band has a nearly square root like
dispersion at the bottom of the band. Fig. 3 shows the comparison of the Al
K-edge spectrum at different temperatures including or excluding explicit
electronic T-dependent effects in the self-energy, namely, using
$\Sigma_{GW}(T=0)$ or $\Sigma_{GW}(T>0)$. When restricting the temperature T
solely to that which is introduced through the density ($\Sigma_{GW}(T=0)$
case), we observe the broadening of the edge and no shift in the edge
position. The shift in the core levels compensates for the shift in the
valence density. On the other hand, for $T$ up to about 1 eV, the finite
temperature self-energy correction is negligible, and the temperature-
independent electron self-energy model is a good approximation. However, as
temperature grows to order $\approx$ 10 eV, the fine structures are smoothed
by the large broadening ($\approx$ 3 eV) associated with shortened electronic
excitation lifetime. The shift in quasiparticle shift due to the finite
temperature self-energy correction is small for the near-edge region, and only
becomes significant between 10 eV and 20 eV above the chemical potential.

As a further illustration of the explicit $T$-dependence of the FT
quasiparticle self-energy, we compute the quasiparticle energy correction
$\Delta_{k}=\varepsilon_{qp}-\varepsilon_{k}$. Here, the quasiparticle energy
$\varepsilon_{qp}$ is the solution to
$\varepsilon_{qp}(k^{\prime})=\varepsilon_{k}+\Sigma_{GW}(k^{\prime},\varepsilon_{qp}(k^{\prime}))$
and $k=\sqrt{2(E-\mu)}$ is the photoelectron wavenumber. We compare the real
part of $\Delta$ in Fig. 4 and imaginary part in Fig. 5 for different self-
energies at the interstitial density: $\Sigma_{GW}(T=0)$ and
$\Sigma_{GW}(T>0)$. Note that the real part of $\Delta$ shows a strong
temperature dependence between 5 eV and 20 eV near the plasmon onset. For the
imaginary part of $\Delta$, the broadening effect becomes important above $T$
= 1 eV.

Figure 3:  K-edge XAS for fcc aluminum ($a=4.05$ Å) using different self-
energies: $T$-independent GW self-energy $\Sigma_{GW}(T=0)$ (solid curves) and
$T$-dependent GW self-energy $\Sigma_{GW}(T>0)$ (dots). Figure 4: The
quasiparticle corrections Re $\Delta_{k}$ for aluminum at temperatures
$T=0,1,2$ and $5$ eV. The calculations used KSDT $v_{xc}(T)$ and different
self-energies: $\Sigma_{GW}(T=0)$ (solid), and $\Sigma_{GW}(T>0)$ (dashed).
Figure 5: The quasiparticle corrections Im $\Delta_{k}$ for aluminum at
temperatures $T=0,1,2$ and $5$ eV. The calculations used KSDT $v_{xc}(T)$ and
different self-energies: $\Sigma_{GW}(T=0)$ (solid), and $\Sigma_{GW}(T>0)$
(dashed).

As a second example, we present results for a noble transition metal (fcc
copper, lattice constant $a=3.61$ Å[54]) for which the $d$-bands are
essentially full. Unlike the K-edge, the L-edge probes the highly-localized
$d$-bands of Cu near the chemical potential. At high temperatures, the pre-
edge peak increases in amplitude due to the decreasing $d$-state occupation
[8, 36]. Fig. 6 shows the $\textnormal{L}_{3,2}$-edge XAS up to $T=5$ eV. The
temperature dependence of $\textnormal{Im}\ \Sigma_{GW}(T)$ results in changes
to pre-edge peaks at $T\gtrapprox 2$ eV. Consequently, the estimation of
temperature based on the pre-edge area method or direct spectrum fitting will
deviate more from the $\Sigma_{GW}(T=0)$ model as temperature increases.

Figure 6: $\textnormal{L}_{3,2}$-edge XAS for Cu ($a=3.61$ Å) finite
electronic temperature $T=0,1,2$ and $5$ eV, where the structure reflects that
of the unfilled $d$-bands. The solid curves denote the $\Sigma_{GW}(T=0)$
self-energy results while the dashes represent the $\Sigma_{GW}(T>0)$ self-
energy results.

##  IV Summary, Conclusions, and Outlook

Our parametrization of the FT GW electron self-energy enables efficient
calculations of XAS at finite $T$, from LT of a few hundred K, up to the WDM
regime with $T$ at least 10 eV. Our strategy uses the QPLDA $G_{0}W_{0}$ level
of refinement for the self-energy, with the RPA dielectric function in
conjunction with the KSDT finite-$T$ LDA exchange-correlation functional.
Specifically, the FT self-energy for a system is approximated using the
uniform electron gas with density equal to that of the local density. This is
a significant simplification, as direct calculations using the exact loss
function of the system for the entire energy range of typical XAS experiments
would be computationally formidable. A finite-temperature SCF procedure for
the XAS calculations is carried out in the complex energy plane in terms of
the FT one-electron Green’s function. The procedure includes the FT exchange-
correlation potential, approximated here by the KSDT parametrization.
Important FT XAS effects include the smearing of the absorption edge and the
presence of peaks in the spectrum below the $T=0$ K Fermi energy. The FT
exchange-correlation potential has only a small effect on XAS at low
temperatures $T\ll T_{F}$ compared to the effect of Fermi smearing. The FT
self-energy is also important for XAS, accounting for both temperature
dependent shifts and final-state broadening. To illustrate its efficacy, the
approach was applied to calculations of XANES for crystalline Al and Cu at
normal density. Above $T>1$ eV, the fine structures experience substantial
broadening in the K-edges, corresponding to a reduction of the quasiparticle
lifetime with increasing $T$.

Going forward, a computationally efficient approximation beyond the uniform
electron gas dielectric function would be to use a many-pole model [56, 57],
which is an extension of the Hedin-Lundqvist single plasmon-pole model [58,
59]. Ab-initio dielectric functions also can be obtained from modern
electronic structure codes. Such a finite-temperature generalization of the
many-pole model is currently under development.

###### Acknowledgements.

The contributions from JJK and JJR are supported by the Theory Institute for
Materials and Energy Spectroscopies (TIMES) at SLAC funded by the U.S. DOE,
Office of Basic Energy Sciences, Division of Materials Sciences and
Engineering under contract EAC02-76SF0051. TST and SBT are supported by DOE
grant DE-SC0002139.

##  Appendix A Model for $\textnormal{Re}\ \Sigma_{GW}$

The real part of the finite temperature GW self-energy, $\textnormal{Re}\
\Sigma_{GW}(T)$, is parametrized using low-order polynomials. The parametric
variables are the Wigner-Seitz radius $r_{s}$, reduced momentum $x=k/k_{F}$
and reduced temperature $t=T/T_{F}$. Below the variable $X$ represents the
array denoted $X=(r_{s},x,t)$.

$\displaystyle\frac{\textnormal{Re}\ \Sigma_{GW}(X)}{E_{F}}$ $\displaystyle=$
$\displaystyle\begin{cases}\alpha_{0}(t)+\alpha_{1}(t)\beta_{1}(r_{s})x\\\
+\alpha_{2}(t)\beta_{2}(r_{s})x^{3/2}\\\
+\alpha_{3}(t)\beta_{3}(r_{s})x^{2}\\\
+\alpha_{4}(t)\beta_{4}(r_{s})x^{5/2},x<\kappa(r_{s},t)\\\ \\\
\tilde{\alpha}_{1}(t)\tilde{\beta}_{1}(t)/x^{-1}\\\
+\tilde{\alpha}_{2}(t)\tilde{\beta}_{2}(t)/x^{-2}\\\
+\tilde{\alpha}_{3}(t)\tilde{\beta}_{3}(t)/x^{-3},x\geq\kappa(r_{s},t)\end{cases}$
(11)

Here $\alpha_{i}(t)=\alpha_{i1}+\alpha_{i2}t+\alpha_{i3}t^{2}$ and
$\tilde{\alpha}_{i}(t)$ has the same form. Similarly
$\beta_{i}(r_{s})=\beta_{i1}r_{s}+\beta_{i1}r_{s}^{3/2}+\beta_{i1}r_{s}^{2}$
and $\tilde{\beta}_{i}(r_{s})$ has the same form. The function $\kappa$ is
defined as:

$\displaystyle\kappa(r_{s},t)$ $\displaystyle=$
$\displaystyle\bigg{(}1+\tanh\bigg{[}a_{1}r_{s}-\alpha_{1}(t)\bigg{]}\bigg{)}$
(12)
$\displaystyle\times\bigg{(}a_{2}\log\big{[}r_{s}\big{]}^{2}+\alpha_{2}(t)\bigg{)}+\alpha_{3}(t)$

where it is fitted to the position of the cusp defined as

cusp $\displaystyle=$
$\displaystyle\begin{cases}\min(1.5,\operatorname*{arg\,min}_{x}\textnormal{Re}\
\Sigma_{GW}),\ 0.2<r_{s}\leq 5\\\
\min(1.5,\operatorname*{arg\,max}_{x}\partial_{x}\textnormal{Re}\
\Sigma_{GW}),\ r_{s}>5\\\ 1,\ r_{s}\leq 0.2\\\ \end{cases}$ (13)

The resulting $\textnormal{Re}\ \Sigma_{GW}$ parametrization is shown in Fig.
7. The absolute mean error for our parametrization is $\approx$ 0.007 $E_{F}$.

Figure 7:  Our parametrization (solid) to the real part of $\Sigma_{GW}(T)$
(gray, dash-dotted) for $r_{s}$ = 0.02, 0.5, 1.0 and 7.0 (blue to red).

##  Appendix B Model for $\textnormal{Im}\ \Sigma_{GW}$

The imaginary part of the finite temperature GW self-energy, $\textnormal{Im}\
\Sigma_{GW}(T)$, is parametrized by

$\displaystyle\frac{\textnormal{Im}\ \Sigma_{GW}(X)}{E_{F}}$ $\displaystyle=$
$\displaystyle\begin{cases}\eta_{1}(r_{s})kt+\eta_{2}(r_{s})kt^{\frac{3}{2}}\\\
+\eta_{3}(r_{s})k^{2}t+\eta_{4}(r_{s})k^{2}t^{\frac{3}{2}}\\\
+\eta_{5}(r_{s})k^{\frac{3}{2}}t+\eta_{6}(r_{s})k^{\frac{3}{2}}t^{\frac{3}{2}}\\\
+\eta_{7}(r_{s})k+\eta_{8}(r_{s})k^{2}\\\ +\eta_{9}(r_{s})k^{\frac{3}{2}}\\\
+\frac{\textnormal{Im}\ \Sigma_{GW}(r_{s},x=1,t)}{E_{F}},x<\lambda(r_{s},t)\\\
\\\
\sigma\sum_{i=1}^{6}\tilde{\eta}_{i}(r_{s})\beta_{i}(t)k^{-i},x\geq\lambda(r_{s},t)\end{cases}$
(14)

where
$\eta_{i}(r_{s}),\tilde{\eta}(r_{s})=\eta_{i1}r_{s}+\eta_{i2}r_{s}^{3/2}+\eta_{i3}r_{s}^{2}$
and $\sigma$ is the standard deviation of the data points used. The function
$\lambda(r_{s},t)$ is given by

$\displaystyle\lambda(r_{s},t)$ $\displaystyle=$
$\displaystyle\bigg{(}1+\tanh\bigg{[}\frac{r_{s}-\frac{1}{2}}{p_{1}}\bigg{]}\bigg{)}$
(15) $\displaystyle\times\bigg{(}p_{2}r_{s}+\eta_{1}(t)\bigg{)}+\eta_{2}(t)$

and Im $\Sigma$ at the Fermi level, $k_{F}$, is parametrized by:

$\displaystyle\frac{\textnormal{Im}\ \Sigma_{GW}(r_{s},x=1,t)}{E_{F}}$
$\displaystyle=$ $\displaystyle\nu_{1}(t)r_{s}^{\frac{1}{2}}+\nu_{2}(t)r_{s}$
(16)
$\displaystyle+\nu_{3}(t)r_{s}^{\frac{3}{2}}+\beta_{4}(t)r_{s}^{\frac{5}{2}}$

where $\nu_{i}(t)=\nu_{i1}t+\nu_{i2}t^{3/2}+\nu_{i3}t^{2}$. The fittings of
Eq. (14), Eq. (15) and Eq. (16) are done for two temperature regions: $t<0.5$
and $t\geq 0.5$, and two density regions: $r_{s}<0.2$ and $r_{s}\geq 0.2$.

The final model is a linear combination of the left and right regions in Eq.
(14). The left region is weighted by the Fermi function
$w=1/\\{1+\exp[15\big{(}k-1.1\lambda(r_{s},t)\big{)}]\\}$ and the right side
by $1-w$. In addition, to prevent spurious negative values at very low $t$, we
clip any negative values to be zero. The resulting parametrization is shown in
Fig. 8. The absolute mean error for our parametrization is $\approx$ 0.006
$E_{F}$.

Figure 8:  Our parametrization (solid) to the imaginary part of
$\Sigma_{GW}(T)$ (gray, dash-dotted) for $r_{s}$ = 0.02, 0.5, 1.0 and 7.0
(blue to red).

## References

  * Peyrusse [2008] O. Peyrusse, Journal of Physics: Condensed Matter 20, 195211 (2008).
  * Recoules and Mazevet [2009] V. Recoules and S. Mazevet, Phys. Rev. B 80, 064110 (2009).
  * Cho  _et al._ [2011] B. I. Cho, K. Engelhorn, A. A. Correa, T. Ogitsu, C. P. Weber, H. J. Lee, J. Feng, P. A. Ni, Y. Ping, A. J. Nelson, D. Prendergast, R. W. Lee, R. W. Falcone, and P. A. Heimann, Phys. Rev. Lett. 106, 167601 (2011).
  * Dorchies  _et al._ [2015] F. Dorchies, V. Recoules, J. Bouchet, C. Fourment, P. M. Leguay, B. I. Cho, K. Engelhorn, M. Nakatsutsumi, C. Ozkan, T. Tschentscher, M. Harmand, S. Toleikis, M. Störmer, E. Galtier, H. J. Lee, B. Nagler, P. A. Heimann, and J. Gaudin, Phys. Rev. B 92, 144201 (2015).
  * Engelhorn  _et al._ [2015] K. Engelhorn, V. Recoules, B. I. Cho, B. Barbrel, S. Mazevet, D. M. Krol, R. W. Falcone, and P. A. Heimann, Phys. Rev. B 91, 214305 (2015).
  * Ogitsu  _et al._ [2018] T. Ogitsu, A. Fernandez-Pañella, S. Hamel, A. A. Correa, D. Prendergast, C. D. Pemmaraju, and Y. Ping, Phys. Rev. B 97, 214203 (2018).
  * Bolis  _et al._ [2019] R. Bolis, J.-A. Hernandez, V. Recoules, M. Guarguaglini, F. Dorchies, N. Jourdain, A. Ravasio, T. Vinci, E. Brambrink, N. Ozaki, J. Bouchet, F. Remus, R. Musella, S. Mazevet, N. J. Hartley, F. Guyot, and A. Benuzzi-Mounaix, Physics of Plasmas 26, 112703 (2019).
  * Jourdain  _et al._ [2020] N. Jourdain, V. Recoules, L. Lecherbourg, P. Renaudin, and F. Dorchies, Phys. Rev. B 101, 125127 (2020).
  * Li  _et al._ [2021] Z. Li, W.-J. L. Li, C. Wang, D. Li, W. Kang, X.-T. He, and P. Zhang, Chinese Physics B 30, 057102 (2021).
  * Faleev  _et al._ [2006] S. V. Faleev, M. van Schilfgaarde, T. Kotani, F. Léonard, and M. P. Desjarlais, Phys. Rev. B 74, 033101 (2006).
  * Hollebon  _et al._ [2019] P. Hollebon, O. Ciricosta, M. P. Desjarlais, C. Cacho, C. Spindloe, E. Springate, I. C. E. Turcu, J. S. Wark, and S. M. Vinko, Phys. Rev. E 100, 043207 (2019).
  * Cho  _et al._ [2016] B. I. Cho, T. Ogitsu, K. Engelhorn, A. A. Correa, Y. Ping, J. W. Lee, L. J. Bae, D. Prendergast, R. W. Falcone, and P. A. Heimann, Scientific Reports 6, 18843 (2016).
  * Denoeud  _et al._ [2014] A. Denoeud, A. Benuzzi-Mounaix, A. Ravasio, F. Dorchies, P. M. Leguay, J. Gaudin, F. Guyot, E. Brambrink, M. Koenig, S. Le Pape, and S. Mazevet, Phys. Rev. Lett. 113, 116404 (2014).
  * Karasiev  _et al._ [2016] V. V. Karasiev, L. Calderín, and S. B. Trickey, Phys. Rev. E 93, 063207 (2016).
  * Karasiev  _et al._ [2018] V. V. Karasiev, J. W. Dufty, and S. Trickey, Phys. Rev. Lett. 120, 076401 (2018).
  * Ichimaru  _et al._ [1987] S. Ichimaru, H. Iyetomi, and S. Tanaka, Physics Reports 149, 91 (1987).
  * Tanaka [2017] S. Tanaka, Contributions to Plasma Physics 57, 126 (2017).
  * Tanaka [2016] S. Tanaka, The Journal of Chemical Physics 145, 214104 (2016).
  * Perrot and Dharma-wardana [2000] F. Perrot and M. W. C. Dharma-wardana, Phys. Rev. B 62, 16536 (2000).
  * Karasiev  _et al._ [2014] V. V. Karasiev, T. Sjostrom, J. Dufty, and S. B. Trickey, Phys. Rev. Lett. 112, 076403 (2014).
  * Groth  _et al._ [2017] S. Groth, T. Dornheim, T. Sjostrom, F. D. Malone, W. M. C. Foulkes, and M. Bonitz, Phys. Rev. Lett. 119, 135001 (2017).
  * Bunău and Calandra [2013a] O. Bunău and M. Calandra, Phys. Rev. B 87, 205105 (2013a).
  * Seah and Dench [1979] M. P. Seah and W. A. Dench, Surface and Interface Analysis 1, 2 (1979).
  * Cushing  _et al._ [2018] S. K. Cushing, M. Zürch, P. M. Kraus, L. M. Carneiro, A. Lee, H.-T. Chang, C. J. Kaplan, and S. R. Leone, Structural Dynamics 5, 054302 (2018).
  * Rehr and Albers [2000] J. J. Rehr and R. C. Albers, Rev. Mod. Phys. 72, 621 (2000).
  * Korringa [1947] J. Korringa, Physica 13, 392 (1947).
  * Kohn and Rostoker [1954] W. Kohn and N. Rostoker, Phys. Rev. 94, 1111 (1954).
  * Dupree [1961] T. H. Dupree, Annals of Physics 15, 63 (1961).
  * Beeby and Edwards [1967] J. L. Beeby and S. F. Edwards, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 302, 113 (1967).
  * Morgan [1966] G. J. Morgan, Proceedings of the Physical Society 89, 365 (1966).
  * Casida [1995] M. Casida, Phys. Rev. A 51, 2005 (1995).
  * Benedict  _et al._ [2002] L. X. Benedict, C. D. Spataru, and S. G. Louie, Phys. Rev. B 66, 085116 (2002).
  * Kas and Rehr [2017] J. J. Kas and J. J. Rehr, Phys. Rev. Lett. 119, 176403 (2017).
  * Allen and Mitrović [1983] P. B. Allen and B. Mitrović (Academic Press, 1983) pp. 1–92.
  * Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
  * Tan  _et al._ [2021] T. S. Tan, J. J. Kas, and J. J. Rehr, Phys. Rev. B 104, 035144 (2021).
  * Martin  _et al._ [2016] R. M. Martin, L. Reining, and D. M. Ceperley,  _Interacting Electrons: Theory and Computational Approaches_ (Cambridge University Press, 2016).
  * Bunău and Calandra [2013b] O. Bunău and M. Calandra, Phys. Rev. B 87, 205105 (2013b).
  * Taillefumier  _et al._ [2002] M. Taillefumier, D. Cabaret, A.-M. Flank, and F. Mauri, Phys. Rev. B 66, 195107 (2002).
  * Gougoussis  _et al._ [2009] C. Gougoussis, M. Calandra, A. P. Seitsonen, and F. Mauri, Phys. Rev. B 80, 075102 (2009).
  * Sham and Kohn [1966] L. J. Sham and W. Kohn, Phys. Rev. 145, 561 (1966).
  * Mustre de Leon  _et al._ [1991] J. Mustre de Leon, J. J. Rehr, S. I. Zabinsky, and R. C. Albers, Phys. Rev. B 44, 4146 (1991).
  * Ankudinov  _et al._ [1998] A. L. Ankudinov, B. Ravel, J. J. Rehr, and S. D. Conradson, Phys. Rev. B 58, 7565 (1998).
  * Mahan [2000] G. Mahan,  _Many-Particle Physics_ (Springer, 2000).
  * Arista and Brandt [1984] N. R. Arista and W. Brandt, Phys. Rev. A 29, 1471 (1984).
  * Kas  _et al._ [2021] J. J. Kas, F. D. Vila, C. D. Pemmaraju, T. S. Tan, and J. J. Rehr, Journal of Synchrotron Radiation 28, 1801 (2021).
  * Rehr  _et al._ [2010] J. J. Rehr, J. J. Kas, F. D. Vila, M. P. Prange, and K. Jorissen, Phys. Chem. Chem. Phys. 12, 5503 (2010).
  * Karasiev  _et al._ [2019] V. V. Karasiev, S. B. Trickey, and J. W. Dufty, Phys. Rev. B 99, 195134 (2019).
  * [49] The FLEUR code (version max-6.0): www.flapw.de.
  * Weinert  _et al._ [1982] M. Weinert, E. Wimmer, and A. J. Freeman, Phys. Rev. B 26, 4571 (1982).
  * Alekseeva  _et al._ [2018] U. Alekseeva, G. Michalicek, D. Wortmann, and S. Blügel, in  _Euro-Par 2018: Parallel Processing_ , edited by M. Aldinucci, L. Padovani, and M. Torquati (Springer International Publishing, Cham, 2018) pp. 735–748.
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Rehr  _et al._ [1994] J. J. Rehr, C. H. Booth, F. Bridges, and S. I. Zabinsky, Phys. Rev. B 49, 12347 (1994).
  * Kittel [2005] C. Kittel,  _Introduction to solid state physics_ , 8th ed. (Wiley, Hoboken, NJ, 2005) p. 20.
  * Kiyono  _et al._ [1978] S. Kiyono, S. Chiba, Y. Hayasi, S. Kato, and S. Mochimaru, Japanese Journal of Applied Physics 17, 212 (1978).
  * Kas  _et al._ [2007] J. J. Kas, A. P. Sorini, M. P. Prange, L. W. Cambell, J. A. Soininen, and J. J. Rehr, Phys. Rev. B 76, 195116 (2007).
  * Kas  _et al._ [2009] J. J. Kas, J. Vinson, N. Trcera, D. Cabaret, E. L. Shirley, and J. J. Rehr, Journal of physics. Conference series 190, 10.1088/1742 (2009).
  * Lundqvist [1967] B. I. Lundqvist, Physik der kondensierten Materie 6, 206 (1967).
  * Hedin and Lundqvist [1971] L. Hedin and B. I. Lundqvist, Journal of Physics C: Solid State Physics 4, 2064 (1971).