# A study of the agreement between binary neutron star ejecta models derived
from numerical relativity simulations

Amelia Henkel  Department of Physics, University of New Hampshire, 9 Library
Way, Durham NH 03824, USA  Francois Foucart  Department of Physics, University
of New Hampshire, 9 Library Way, Durham NH 03824, USA  Geert Raaijmakers
GRAPPA Institute of High-Energy Physics, University of Amsterdam, Science Park
904,
1098 XH Amsterdam, Netherlands  Samaya Nissanke  GRAPPA Institute of High-
Energy Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, Netherlands

###### Abstract

Neutron star mergers have recently become a tool to study extreme gravity,
nucleosynthesis, and the chemical composition of the Universe. To date, there
has been one joint gravitational and electromagnetic observation of a binary
neutron star merger, GW170817, as well as a solely gravitational observation,
GW190425. In order to accurately identify and interpret electromagnetic
signals of neutron star mergers, better models of the matter outflows
generated by these mergers are required. We compare a series of ejecta models
to see where they provide strong constraints on the amount of ejected mass
expected from a system, and where systematic uncertainties in current models
prevent us from reliably extracting information from observed events. We also
examine 2396 neutron star equations of state compatible with GW170817 to see
whether a given ejecta mass could be reasonably produced with a neutron star
of said equation of state, and whether different ejecta models provide
consistent predictions. We find that the difference between models is often
comparable to or larger than the error generally assumed for these models,
implying better constraints on the models are needed. We also note that the
extrapolation of outflow models outside of their calibration window, while
commonly needed to analyze gravitational wave events, is extremely unreliable
and occasionally leads to completely unphysical results.

##  I Introduction

The gravitational observation of binary neutron star merger GW170817 and its
corresponding kilonova signal AT2017gfo ushered the world into the era of
multi-messenger astronomy (Abbott et al. 2017 [1]). Binary neutron star (BNS)
mergers in general are interesting for the wealth of information encoded in
their signal. The electromagnetic signal of BNS mergers can elucidate several
source parameters, such as the mass of the outflows from the system, its
composition, morphology, and potentially information about the fate of the
system. Their gravitational wave signal can provide information about the
masses of the component objects, their mass ratio, spins, and the fate of the
remnant system. In general, binary neutron star mergers contain rich physics,
providing insight into the strong gravity regime; nucleosynthesis and the
origin of heavy elements; the dynamics and formation mechanisms of gamma ray
bursts (GRBs); and the neutron star equation of state (EOS).

The ejecta and observed kilonova signal is highly dependent on the mass,
velocity, morphology, and composition of the outflows (Wollaeger et al. 2021
[2], Barnes and Kasen 2013 [3]). We often model the merger ejecta as
consisting of two components: the dynamical ejecta, ejected within
milliseconds of the merger, and the disk wind ejecta, produced up to a few
seconds after the merger. The dynamical ejecta is typically faster, and may
include a neutron-rich component associated with the tidal disruption of a
neutron star and a less neutron-rich component due to shocks and oscillations
in the forming neutron star remnant. The wind ejecta is typically slower,
though both its velocities and composition may be significantly impacted by
physical processes that are not fully included in many existing post-merger
simulations (magnetohydrodynamics, neutrino transport) [4, 5, 6].

Given the limited number of neutron star merger simulations performed so far,
candidate electromagnetic signals to neutron star mergers are typically
analyzed using approximate analytical model fitted to the result of numerical
simulations. To date, a number of models have been developed to constrain the
dynamics and outflows of binary neutron star mergers (e.g. [7, 8, 9, 10, 2,
11]). Recent work by Nedora et al. 2021 [11] clearly demonstrated that
analytical fits attempting to predict the properties of neutron star merger
ejecta can be heavily influenced by the level of microphysics included in the
simulations used to derive that fit. In this manuscript, we attempt to answer
a slightly different question, namely how robust inferences made about the
properties of neutron stars are to the choice of fitting formula. This remains
an important open question because some models do not include reliable error
bars for their predictions, while models that do include errors bars only
cover a limited region of parameter space and may lead to significantly larger
systematic errors when used outside of their intended region of validity –
something that should ideally be avoided, yet is common practice when
analyzing existing EM signals from neutron star mergers due to the lack of
models covering the entire parameter space of interest.

We structure our paper in the following way: first, we introduce the ejecta
models considered in this study. We then present a comparison of these models,
and show where they currently agree and disagree, under the simplifying
assumption that both neutron stars have the same radius. We later drop the
constant radius assumption and apply the ejecta models to a series of
equations of state consistent with gravitational wave observations of GW170817
to see if current observations can elucidate the true dense matter EOS.
Finally, we conclude with a discussion of our results, drawing attention to
the cases most similar to BNS events GW170817 and GW190425.

##  II Ejecta models

Source | Model | Simulations | Inputs | Calibrated Q | Uncertainty
---|---|---|---|---|---
DU [8] | Dyn. Ej. | 172 sims from [12, 13, 14, 15, 16, 17] | $M_{1},M_{2},C_{1},C_{2}$ | 0.48-1.0 | $0.004M_{\odot}$
KF [10] | Dyn. Ej. | 200 sims from [8, 18] | $M_{1},M_{2},C_{1},C_{2}$ | 0.48-1.0 | $0.004M_{\odot}$
NAL [11] | Dyn. Ej. | 324 sims from [11, 19, 16, 20, 9, 21, 22, 23] | $\tilde{\Lambda},\tilde{Q}$ | 0.54-1.0 | $0.5M_{ej}+5\times 10^{-5}M_{\odot}$
KF [10] | Disk | 57 sims from [18, 9] | $M_{low},C_{low}$ | 0.775-1.0 | $0.5M_{disk}+5\times 10^{-4}M_{\odot}$
DAL [24] | Disk | 73 sims from [18, 9, 17, 25] | $Q,M_{tot},M_{th}$ | 0.571-1.0 | Not quoted
RAL [9] | Disk | 59 sims performed by RAL | $\tilde{\Lambda}$ | 0.85-1.0 | $0.5M_{disk}+5\times 10^{-4}M_{\odot}$
NAL [11] | Disk | 119 sims from [11, 19, 16, 20, 9, 21, 22, 23] | $\tilde{\Lambda},\tilde{Q}$ | 0.54-1.0 | $0.5M_{disk}+(5\times 10^{-4})M_{\odot}$
Table 1: Information about the dynamical ejecta and disk mass models
considered in this work. From left to right, the columns indicate (i) the name
of the model; (ii) type of outflow; (iii) information about the simulations to
which the models are calibrated; (iv) the models’ input parameters; (v) the
range of mass ratios to which the models are calibrated; and (vi) the
uncertainty associated with each model, when quoted.

###  II.1 Dynamical Ejecta Models

First, we visit three fitting formulae for the amount of dynamical ejecta
emitted from neutron star mergers, taken from Kruger and Foucart 2020 [10]
(hereafter KF), Dietrich and Ujevic [8] (hereafter DU), and Nedora et al. [11]
(hereafter NAL). We note that these formulae do not distinguish between the
cold, neutron-rich tidal ejecta and the hot, less neutron-rich ejecta produced
by the collision of the neutron star cores, and do not tell us anything about
the geometry of the ejecta. We do not investigate here the impact of these
issues on kilonovae light curves, but it is worth noting that even if these
formulae were exact, they would not be sufficient to predict the kilonova
signals associated with the dynamical ejecta. NAL introduced some fitting
formulae providing partial information about the composition and geometry of
the outflows, but they point out that the accuracy of these formulae is
significantly limited by the small number of numerical relativity simulations
with sufficiently advanced microphysics that reported the observables fitted
in their work.

DU present a phenomenological fit of the mass, energy, and velocities of
dynamical ejecta derived from a large series of numerical relativity
simulations (N=172). They use simulations presented in Hotokezaka et al. 2013
[12], Bauswein et al. 2015 [13], Dietrich et al. 2015 [14], Lehner et al. 2016
[15], Sekiguchi et al. 2016 [16], and Dietrich et al. 2017b [17]. Their data
set combines results from general relativistic grid structured codes with
results employing a smoothed particle hydrodynamics (SPH) code under the
conformal flatness approximation. It includes simplified EOSs and tabulated
EOSs, as well as simulations with and without neutrino treatment. As this
study only considers dynamical ejecta, and no wind ejecta, Dietrich and Ujevic
note that their estimates can provide a lower bound for the luminosity of EM
observables.

Their fit for the dynamical ejecta mass is as follows:

$\frac{M_{ej}^{fit}}{10^{-3}M_{\odot}}=\left[a\left(\frac{M_{2}}{M_{1}}\right)^{1/3}\left(\frac{1-2C_{1}}{C_{1}}\right)\right.\\\
\left.+b\left(\frac{M_{2}}{M_{1}}\right)^{n}+c\left(1-\frac{M_{1}}{M_{1}^{*}}\right)\right]+1\leftrightarrow
2+d$ (1)

where $a=-1.35695$, $b=6.11252$, $c=-49.43355$, $d=16.1144$, and $n=-2.5484$;
$M^{*}$ is the baryonic mass, M is the gravitational (ADM) mass in isolation,
$C=GM/Rc^{2}$ is the compactness and $R$ is the radius of the neutron star.
KF developed formulae for the dynamical ejecta of BNS and black hole-neutron
star (BHNS) systems, as well as the disk mass for BNSs. Their dynamical ejecta
fit for BHNS systems is outside the scope of this work, as we only consider
BNSs; their disk mass ejecta fit is described in greater detail in Section
2.2. KF’s fit is similar to DU’s, but does not depend on baryonic mass. In
their fitting formulae, they aim for analytical simplicity as well as
physically reasonable extrapolation towards high-compactness stars. They also
derive their fits from a slightly larger region of parameter space; namely,
they base their work off of NR simulations from Radice et al. 2018 [9] and
Kiuchi et al. 2019 [18], in addition to those used by DU. These additional
simulations include, among other things, more asymmetric binaries than in the
original data set.

Their dynamical ejecta fit for binary neutron stars is given as

$\frac{M_{dyn}}{10^{-3}M_{\odot}}=\left(\frac{a}{C_{1}}+b\frac{M_{2}^{n}}{M_{1}^{n}}+cC_{1}\right)M_{1}+(1\leftrightarrow
2)$ (2)

Here, a=-9.3335, b=114.17, c=-337.56, and n=1.5465.

NAL presents fitting formulae for the dynamical ejecta of BNS mergers based on
a polynomial in mass ratio $\tilde{Q}$ and reduced tidal deformability
$\tilde{\Lambda}$, defined as

$\tilde{\Lambda}=\frac{16}{13}\frac{(M_{A}+12M_{B})M_{A}^{4}\Lambda_{A}}{M^{5}}+(A\leftrightarrow
B).$ (3)

In this convention, $\tilde{Q}=M_{A}/M_{B}\geq 1$ and $M=M_{A}+M_{B}$. They
obtain their fitting formula from a suite of recent NR simulations (available
at Nedora et al. 2020 [26]) with varying levels of realism in their treatment
of the composition of the star and neutrino transport. While they present
fitting formulae with two distinct polynomial orders and including or
excluding simulations with better/worse microphysics, that they then apply to
a multitude of binary parameters, we only consider here their “recommended”
fitting formula applied to the dynamical ejecta mass and disk mass. The
fitting formula is as follows:

$P_{2}^{2}(\tilde{Q},\tilde{\Lambda})=b_{0}+b_{1}\tilde{Q}+b_{2}\tilde{\Lambda}+b_{3}q^{2}+b_{4}\tilde{Q}\tilde{\Lambda}+b_{5}\tilde{\Lambda}^{2}.$
(4)

For the dynamical ejecta, they provide best-fit parameters $b_{0}=-1.32$,
$b_{1}=-0.382$, $b_{2}=-4.47\times 10^{-3}$, $b_{3}=-0.339$, $b_{4}=3.21\times
10^{-3}$, and $b_{5}=4.31\times 10^{-7}$.

###  II.2 Disk Mass Models

Fitting formulae for the mass remaining in an accretion disk around the
remnant black hole (or neutron star) after a BNS merger have also been
provided in multiple works. These formulae are typically fit to a lower number
of simulations than dynamical ejecta formulae, as not all numerical relativity
simulations report remnant disk masses. Additionally, the definition of the
disk mass itself varies between different studies, complicating comparisons
between different sets of simulations. Most notably, reported disk masses are
sensitive to the time at which they are measured, and to the definition of the
boundary between the “disk” and “neutron star” in neutron star-disk systems.
They are however crucial to kilonovae modeling, as most BNS mergers likely
eject more mass through disk winds than dynamical ejecta.

KF provides a fitting formula for the disk mass ejected from a binary neutron
star merger, based on the subset of simulations used for the dynamical ejecta
model that provide disk mass information. The disk mass fit for BNSs they
arrive at is as follows:

$M_{disk}=M_{1}\textrm{max}\\{5\times 10^{-4},(aC_{1}+c)^{d}\\}$ (5)

with $C_{1}$ the compactness of the lower mass neutron star, and best-fit
coefficients a=-8.1324, c=1.4820, and d=1.7784. Their fit has an associated
uncertainty of

$\Delta M_{disk}=0.5M_{disk}+5\times 10^{-4}M_{\odot}$ (6)

We also consider the disk mass fitting formula provided by Dietrich et al.
2020 [24], hereafter DAL. Their fit is derived from a series of 73 numerical
relativity simulations performed by Kiuchi et al. [18], Radice et al. [9],
Dietrich et al. [17], and Hotokezaka et al. [25]. Their model takes a similar
form to their previous results in Coughlin et al. 2019 [27], but improves upon
it by including mass-ratio dependent fitting parameters. They arrive at the
following disk mass fit:

$log_{10}(M_{disk})=max\left(-3,a\right.\\\
\left.\left(1+\beta\textrm{tanh}\left(\frac{c-M/M_{threshold}}{d}\right)\right)\right)$
(7)

where $a$ and $b$ are given as

$a=a_{0}+\delta a\xi$ $b=b_{0}+\delta b\xi$

and $\xi$ can be expressed as

$\xi=\frac{1}{2}\textrm{tanh}\left(\beta(Q-Q_{trans})\right).$ (8)

Here, the best-fit coefficients are $a_{0}$=-1.581, $\delta a$=-2.439,
$b_{0}$=-0.538, $\delta b$=-0.406, c=0.953, d=0.0417, $\beta$=3.910, and
$Q_{trans}$=0.900.

Additionally, we consider the model provided by NAL. Their polynomial fit for
the dynamical ejecta (Equation 4) is also used for the disk mass. For the disk
mass, they have that $b_{0}=-1.85$, $b_{1}=2.59$, $b_{2}=7.07\times 10^{-4}$,
$b_{3}=-0.733$, $b_{4}=-8.08\times 10^{-4}$, and $b_{5}=2.75\times 10^{-7}$.

Lastly, we employ the results from Radice et al. 2018 [9], hereafter RAL. They
study the mass ejection, nucleosynthetic yields, and electromagnetic
counterparts of BNS mergers based on 59 high-resolution numerical relativity
simulations. They employ four nuclear-theory based equations of state (SFHo,
BHB$\Lambda\phi$, DD2, and LS220) to constrain current estimates on merger
outflows and ejecta properties. Their models are calibrated to a mass ratio
range of 0.85-1.0. Their data is fit by the following expression:

$\frac{M_{disk}}{M_{\odot}}=max\left\\{10^{-3},\alpha+\beta\textrm{tanh}\left(\frac{\tilde{\Lambda}-\gamma}{\delta}\right)\right\\}$
(9)

They also quote the uncertainty in the disk mass as

$\Delta M_{disk}=0.5M_{disk}+5\times 10^{-4}M_{\odot}$ (10)

##  III Methods

###  III.1 Comparing Dynamical Ejecta Models

First, we compare dynamical ejecta models from DU, KF, and NAL. We do this to
see for what radii and mass ratios the models agree, and to measure their
level of agreement over the parameter space. For each figure, we fix the chirp
mass, a quantity that is in practice obtained from GW observations. We then
generate a range of 100 mass ratios $Q=M_{1}/M_{2}$ from 0.5 to 1.0 (where
$M_{1}<M_{2}$), which allows us to calculate the individual masses $M_{1}$ and
$M_{2}$ with the following relation:

$M_{2}=M_{c}\frac{(Q+1)^{1/5}}{Q^{3/5}}.$ (11)

In practice, $M_{c}$ will typically be fairly well measured, but $Q$ will come
with significant uncertainties. Additionally, measurements of $Q$ are
partially degenerate with measurements of the spins of the neutron stars. In
this manuscript, we focus on uncertainties due solely to the choice of the
ejecta model, and thus assume fixed values of $M_{c}$ and $Q$; yet in actual
observations the uncertainty in $Q$ would certainly have to be taken into
account.

We then generate 100 values of the neutron star radius from 10km to 14km,
working under the assumption that both neutron stars have the same radius;
$R_{1}=R_{2}$. The use of a constant radius is clearly an approximation. As
evidenced by Figures 4 and 5, however, a broad range of masses correspond to a
very small range of radii for most equations of state consistent with existing
nuclear physics and astrophysical constraints – and visualizing models in the
2D space of mass ratios and radii is easier than in the infinite dimensional
space of mass ratio and equations of state. We will lift this simplifying
assumption later when considering predictions for specific equations of state.

We select chirp masses $M_{c}=1.44M_{\odot}$, the value associated with
GW190425 [28], and $M_{c}=1.186M_{\odot}$, the value associated with GW170817
[1]. After calculating the component masses, the dynamical ejecta mass for
each combination of R and Q is calculated according to Eqns. 1, 2, and 4. We
then plot the dynamical ejecta from each model as color gradients on a radius
vs. mass ratio grid, as seen in Figure 1.

It is evident from Figure 1 that for both the DU and KF dynamical ejecta
models, a higher ejecta mass is predicted for more unequal mass mergers. This
is to be expected, as unequal mass mergers are associated with more tidal
deformation. Additionally, KF’s model is more sensitive to radius for both
cases. One objective of their model was to account for the expected lack of
ejecta from highly compact stars near the equal mass limit, a phenomena that
is physically understood but not automatically reflected by fitting formulae
due to the lack of numerical simulations in that regime. This is particularly
visible in the higher chirp mass case, where a larger fraction of the systems
are expected to collapse to a black hole. The NAL model shows very different
qualitative behavior, especially outside of the regions where most numerical
relativity simulations used to calibrate these models are found (near equal
mass, non-collapsing systems). We believe part of this discrepancy can be
attributed to the fact that NAL’s quadratic ejecta model takes a different
functional form than KF and DU. Additionally, it is a function of $\tilde{Q}$
and $\tilde{\Lambda}$, whereas the other two models are functions of the
component masses and compactnesses. It is important to note that while the KF
model does recover what we expect to be the correct limit for small radii,
there is no evidence that it performs better than the other models at more
moderate radii. For example, Nedora et al. 2021 showed that the NAL and KF
models have, despite their very visible qualitative differences, comparable
fitting residuals when fitted to their chosen datasets, while Camilletti et
al. 2022 [29] showed that existing dynamical ejecta models compare poorly to
numerical simulations for GW190425-like system (with the KF model performing
especially poorly in that regime). There is also no particular reason to
believe that one model is more accurate when extrapolating to very asymmetric
systems – in fact, all models predict unrealistically large ejecta mass for
asymmetric systems, and it is more likely than not that they are all very
inaccurate in that regime (the models predict masses well above
$0.06M_{\odot}$, while no simulation used to calibrate the model has seen
ejected masses above $0.06M_{\odot}$, and most simulations find dynamical
ejecta masses $M_{\rm ej}\lesssim 0.01M_{\odot}$). Using these models for
asymmetric binaries can thus be particularly dangerous if no corrections are
applied, especially when assessing the impact of mergers on r-process
nucleosynthesis (see e.g. [30]).

Figure 1: Dynamical ejecta mass predicted by Dietrich and Ujevic’s model
(left); Kruger and Foucart’s model (center); and Nedora et al.’s model
(right). The contours are in units of solar masses; the chirp mass is
$M_{c}=1.186M_{\odot}$ for the top panel and $M_{c}=1.44M_{\odot}$ for the
lower panel. We truncate values above $0.022M_{\odot}$, to make the more
reliable regions of the parameter space easier to study. At high mass ratio,
the extrapolation of fitting formulae outside of their range of applicability
leads to extremely large predicted masses that have never been observed so far
in simulations $\gtrsim 0.06M_{\odot}$).

We note that while a mass ratio $Q\sim 0.5$ is not physical for the low chirp
mass case (the lower mass star would have $M\sim M_{\odot}$), it is a
perfectly valid possibility for the high chirp mass case (lower mass of $\sim
1.2M_{\odot}$).

###  III.2 Disk Mass Model Comparison

Next, we perform the same procedure with the aforementioned disk mass models:
KF (from [10]), DAL (from [24]), RAL (from [9]), and NAL (from [11]). This was
done to determine if any constraints about a BNS EM observation can be made
given some level of agreement between said models.

Figure 2: Disk mass ejecta calculated with each of the four disk mass models
considered in this study. The figures in the top panel have a chirp mass of
$M_{c}=1.188M_{\odot}$; the figures in the bottom panel have
$M_{c}=1.448M_{\odot}$. We truncate the ejecta values for
$M_{disk}>0.499M_{\odot}$ because the values quickly become unphysical outside
of their domains.

We can immediately observe that some predicted disk masses become unphysical
for increasingly unequal mass cases, similar to the dynamical ejecta. For
example, the NAL model predicts upwards of 16$M_{\odot}$ of disk ejecta for
low compactness (larger R, smaller Q) when $M_{c}=1.188M_{\odot}$. For ease of
visualization of the results, we truncate any values in excess of
0.5$M_{\odot}$ in our figures.
All the disk mass models considered in this work predict little to no disk for
nearly equal mass systems; this is to be expected, as the most efficient way
to form massive accretion disks in BNS mergers is through the tidal disruption
of a lower mass star by its more massive companion. That less compact (larger
R) stars produce more massive disks is well-captured in the DAL, KF, and RAL
models. NAL has a more complex behavior as a function of $R$, especially for
high mass systems, due to the chosen functional form for the model (nearly all
models from NAL are quadratic with a saddle point close to or within the
fitting region).

While BNSs with $Q\sim 0.5$ are highly unlikely given the small range of
possible NS masses, we deliberately probe a large parameter space to ascertain
the models’ behavior at extreme – and less well-studied – regions of the
parameter space.

As for the dynamical ejecta models, the disk mass models show particularly
significant disagreement for GW190425-like systems. Camilletti et al. [29]
showed that for GW190425-like systems, the KF model is accurate (within its
relatively large error bars) in the parts of the high-mass parameter space
covered by their simulations (for $Q\gtrsim 0.65$), though we caution again
about inferring too much from that result about the quality of the model in
other regions of the BNS parameter space.

###  III.3 Equation of State

The previous sections assumed constant-radius equations of state; an
assumption convenient for visualizing the results but that does not exactly
match realistic equations of state. For a more realistic treatment of the
nuclear equation of state, we consider a suite of 2396 neutron star equations
of state (EOSs). The EOSs employed in this study are a series of best-fit
“spectral” equations of state calibrated from GW170817 data. More information
about the EOSs included in this study can be found in Abbott et al. 2018 [31],
while the theoretical framework on which these equations of state are based
can be found in Lindblom 2010 [32]. Each EOS in the publicly available master
file contains best-fit values for the spectral EOS parameters $\gamma_{0}$,
$\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$, as well as the pressure at a
reference density. Information about the mass, radius, and tidal deformability
for each EOS can then be obtained from this data (see Raaijmakers et al. 2020
[33]).
The main practical use of the analytical models discussed here is to attempt
to determine whether the binary parameters of a given system - or a series of
binary parameters from multiple events - can be used to reliably rule out
certain equations of state, and thus constrain the true neutron star EOS. In
this manuscript, we test how robust these inferences are to the choice of
model assuming extremely simplified constraints from potential observations.

First, we define a fiducial system’s ejecta mass, chirp mass and mass ratio.
In practice, the chirp mass and mass ratio are only known up to the
potentially large uncertainties in gravitational wave parameter estimation,
while the derivation of an ejecta mass from e.g. a kilonova signal is a non-
trivial process that depends on the morphology, composition, and velocity of
the ejecta as well as nuclear physics and radiation transport uncertainties.
In order to directly study uncertainties due to the current modeling of the
ejecta mass, however, we ignore these important complications for now, and do
not consider any uncertainty in the fiducial ejecta mass, chirp mass, or mass
ratio. We define the fiducial ejecta mass as the total amount of dynamical
ejecta and unbound disk ejecta; in practice, this quantity would have to be
inferred from kilonova observations (see e.g. Raaijmakers et al. 2021 [33] and
Coughlin et al. 2019 [27]). Defining $M_{c}$ and $Q$ sets the component masses
of our system $M_{1}$ and $M_{2}$. In the plots discussed in this section,
$M_{c}$ and $Q$ are fixed in each subplot, and thus each subplot corresponds
to some constant $M_{1}$ and $M_{2}$.

Second, we obtain the radii $R_{1}$ and $R_{2}$ and tidal deformabilities
$\Lambda_{1}$ and $\Lambda_{2}$ corresponding to $M_{1}$ and $M_{2}$
respectively, for each equation of state. From these values, we can calculate
the compactness $C$ of each neutron star as well as the reduced tidal
deformability $\tilde{\Lambda}$ of the binary. These quantities all serve as
inputs for the ejecta models considered in this work.

Then we calculate the dynamical ejecta as a function of the aforementioned
inputs for each equation of state. Each subplot considers either the model
presented by KF, DU, or NAL; the dynamical ejecta model considered for each
subplot is indicated in its title.

Next, we calculate the disk ejecta, also functions of the aforementioned
quantities, for each equation of state, assuming that a fraction $f_{\rm
disk}$ of the accretion disk is ejected as disk wind in the seconds following
the merger. Given uncertainties $\sigma_{dyn}$ and $\sigma_{disk}$ in
$M_{dyn}$ and $M_{disk}$ and assuming $f_{\rm disk}\in[f_{\rm low},f_{\rm
high}]$, we calculate a “window” of possible values for the total ejected
mass:

$M_{min}=(M_{dyn}-\sigma_{dyn})+\textrm{f}_{low}(M_{disk}-\sigma_{disk})$
$M_{max}=(M_{dyn}+\sigma_{dyn})+\textrm{f}_{high}(M_{disk}+\sigma_{disk})$

Doing so provides a range of ejecta values that, according to the model, is
compatible with the input parameters $M_{c}$, $Q$ and the chosen equation of
state. The fiducial ejecta mass is then compared to each calculated ejecta
mass window to see if the fiducial mass falls within the calculated range. For
the KF and DU figures, we present a comparison of KF, DAL, and RAL’s disk
ejecta models. For the NAL figures, we overlay the performance results of KF,
DAL, RAL, and NAL’s disk ejecta models. This is done to more clearly draw a
comparison between what ejecta models of different functional forms predict.
Additionally, we hold fixed the fractional range of unbound disk mass, with
$f_{low}=0.1$ and $f_{high}=0.4$. We initially considered two windows,
$10\%-40\%$ and $20\%-25\%$, to check whether better constraints on the
fraction of unbound disk mass impacts our results, but found that doing so
does not provide any additional information at this point – the differences
between models and fit uncertainties have a stronger impact on our results.

If for a given equation of state the fiducial ejecta mass falls within the
window of acceptable masses for all models under consideration, then that
equation of state is assumed to be compatible with the fiducial observation.
If it falls outside of that window for all models, the equation of state is
incompatible with that observation. Finally, if the fiducial ejecta mass falls
within the allowed region for some models but not other, no robust inference
can be made about that equation of state – and inferences made without
considering modeling uncertainties may be inaccurate. For each subplot, we
overlay the mass-radius curve of each of the 2396 equations of state onto one
plot, color coded according to whether the fiducial binary parameters agree
with none (black), some (red), or all (green) of the disk mass ejecta models.
We demonstrate this process in Figure 3.

Our results are presented in Figures 4 and 5. We consider the same chirp
masses as before, $1.186M_{\odot}$ and $1.44M_{\odot}$, and the two mass
ratios $Q=0.9$ and $Q=0.72$. We use $Q=0.72$ to demonstrate highly unequal
mass mergers; any mass ratio more extreme produces values of $M_{1}$ and
$M_{2}$ which fall outside the range allowed by some of the EOSs. We ignore
mass ratios $Q\lesssim 0.7$ as we have already seen that the existing models
are unreliable in that regime. Finally, we initialize fiducial ejecta masses
of $0.05M_{\odot}$, which is similar to GW170817, and $0.005M_{\odot}$, to
study which systems are compatible with effectively negligible mass ejection
(such as, possibly, GW190425).

Figure 3: A demonstration for how the subsequent figures in this paper are
produced. For the suite of EOSs considered in this work, we determine whether
each one can reproduce ejecta in agreement with a pre-selected fiducial value,
and color-code their mass-radius relations accordingly. This is done by, for a
given set of input parameters ($M_{c}$ and $Q$), calculating the dynamical
ejecta according to one model, then calculating the disk ejecta according to
multiple models for comparison, then overlaying where the results agree and
disagree. The larger figure shows an overlay of the results from (NAL dyn.
ej.) + (disk ej.) for the four models considered in this work. The four
smaller subplots show the results for each of the disk models.
Figure 4: Mass-radius relationships for 2396 equations of state from [31]. The
figures above employ one of three dynamical ejecta models, KF, DU, or NAL; one
of two fiducial ejecta masses, 0.05 and 0.005; and one of two chirp masses,
$1.186M_{\odot}$ (in agreement with observations of GW170817) and
$1.44M_{\odot}$ (in agreement with GW190425). The KF and DU dyn. ej. figures
overlay a comparison of KF, DAL, and RAL disk ejecta, whereas the NAL dyn. ej.
figures overlay a comparison of KF, DAL, RAL, and NAL disk ejecta. All figures
here correspond to systems with a fixed mass ratio Q=0.9.
Figure 5: The same scenarios as Figure 4, but with mass ratio $Q=0.72$.

We present all possible configurations of the input parameters in Figures 4
and 5, with each row corresponding to the same $M_{c}$, $M_{t}$ and $Q$ for
ease of comparison. Looking first at Figure 4, we can see that for
$M_{c}=1.448M_{\odot}$, the KF and DU models consistently agree that for all
but the stiffest equations of state, it is not feasible for such a system to
produce significant ejecta. For such a system, NAL predicts the opposite:
significant ejecta for all except the stiffest EOSs. As can be seen in Figure
2, DAL, KF, and RAL predict almost no ejecta to be produced for mass ratios
above $Q\sim 0.75$, while NAL’s disk model predicts the presence of disk
ejecta up to around $Q\sim 0.9$. This explains why the NAL figures appear red
for columns 3 and 4 of Figure 4, while KF and DU are both green (for
$M_{t}=0.005M_{\odot}$) and both black (for $M_{t}=0.05M_{\odot}$); NAL’s disk
model is in direct disagreement here with the other models.

However, the models tend to agree in other parts of the parameter space. This
is shown in the second row of Figure 4, which is also the case most similar to
GW170817 ($M_{c}=1.186M_{\odot}$, significant ejecta). Here, only the least
compact EOSs are consistently allowed. We do observe some disagreement for the
less massive EOSs; the KF and DAL subplots are able to consistently rule out
approximately $30\%$ of the least stiff ones, but this not possible for the
NAL figure. This is likely because NAL’s disk mass model produces relatively
more ejecta for NSs with smaller radii compared to KF and DU.

We see that the trend for $M_{c}=1.186,Q=0.9,M_{t}=0.05M_{\odot}$ is inverted
when we consider instead $M_{t}=0.005M_{\odot}$ of ejecta. Here, only the most
compact EOSs are able to predict such little mass. Each subplot disallows a
varying amount of EOSs: KF, DAL, and NAL each reject approximately $35\%$,
$25\%$, and $7\%$ of the least compact EOSs. Because each subplot shows this
least compact $\sim 7\%$ as being disallowed, we can say that these EOSs would
be reliably ruled out by such an observation.

We now turn the focus of our discussion to Figure 5, which corresponds to
$Q=0.72$. Looking first to the top row which corresponds to
$M_{c}=1.186M_{\odot}$ and a fiducial mass $M_{t}=0.005M_{\odot}$, we can see
that KF and DAL subplots agree that such an unequal system cannot produce
negligible ejecta. The narrative shifts for the NAL subplot, where NAL’s low
dynamical ejecta appears to be the source of disagreement. We can see from
Figure 1 that while DAL and KF predict significant dynamical ejecta for nearly
all radii in this regime, NAL predicts almost none. The second row of Figure
5, which corresponds also to $M_{c}=1.186M_{\odot}$ but now a fiducial mass
$M_{t}=0.05M_{\odot}$, demonstrates surprisingly consistent results. (A close
inspection shows that the NAL subplot has slightly less green-colored EOSs;
the margin at which all models begin to agree occurs slightly to the right of
where the margin exists for KF and DU.) While this case does not robustly rule
out certain EOSs, it does consistently identify viable ones – here, the least
compact EOSs are allowed. The fourth row of Figure 5, with
$M_{c}=1.44M_{\odot}$ and $M_{t}=0.05M_{\odot}$ of ejecta, shows significant
disagreement between ejecta models. While the KF subplot disallows the most
compact EOSs, DU and NAL do not arrive at this result. That each subplot has
around half the EOSs colored red in the same region is not informative; if
anything, this demonstrates the models are consistently unable to make robust
predictions in this region of parameter space. Similarly, we would be unable
to draw robust conclusions from an observed system with $M_{c}=1.44M_{\odot}$
and $M_{t}=0.005M_{\odot}$. When using DAL’s dynamical ejecta model, all EOSs
are disallowed; when using KF’s dynamical ejecta model, we see the most
compact EOSs allowed and least compact are consistently disallowed; and when
using NAL’s dynamical ejecta model, the models disagree for the entire
parameter space. While the results here disagree with each other, we can see
that they are consistent with what is predicted in Figure 1. DAL’s dynamical
ejecta model produces more mass at $Q=0.72$ for all radii, disqualifying all
EOSs in the subplot. We can also see from Figure 2 that NAL predicts, for most
radii, a vastly larger disk than the other three disk models.

##  IV Discussion

We have considered a series of dynamical ejecta and disk mass models for
binary neutron star mergers and estimated their predictive power by measuring
agreement between models for a given system. We find that regions where all
existing models are consistent are quite rare, and that we should thus be
cautious about making inferences about the parameters of BNS systems using
just one of these models – even when accounting for the models’ stated
uncertainties. We recommend to at least compare results with different ejecta
models in order to get a sense of modeling uncertainties for a given system.

In more detail, we find that for the dynamical ejecta models, DU and KF
demonstrate a stronger dependence on mass ratio than is seen for NAL. That
more unequal mass systems produce relatively more ejecta compared to equal
mass systems is well-captured in DU and KF’s models, but not by NAL. We note
that this is likely explained by NAL’s model being a function of different
parameters compared to DU and KF ($Q$ and $\tilde{\Lambda}$ vs. $M_{1,2}$ and
$C_{1,2}$), and that the model takes a different functional form compared to
the two other models. Of the three models, KF demonstrates the strongest
dependence on radius, which is an intentional feature meant to capture the
different behavior of BNS systems of varying compactness. All of the models
begin to diverge to unrealistically high values as we probe increasingly
unequal mass systems, which is due to the models being used outside of the
parameter space to which they were calibrated. For this reason, it can often
be problematic to extrapolate ejecta fits. We also wish to highlight that the
differences between what the models predict is often of the order of the
actual predicted ejecta mass.

We perform the same study on the disk mass models considered in this work. The
models consistently predict low/no disk masses for compact (small R), equal-
mass systems, which agrees with our current understanding that such mergers
result in prompt collapse to a black hole. However, different models have
different thresholds for when ejecta is or is not produced. All disk models
with the exception of NAL produce significant ejecta for increasingly unequal
mass systems; for $Q\sim 0.55-0.6$, the model predicts negligible disk mass.
As was the case for the dynamical ejecta, we see that the differences between
the models are comparable in magnitude to the models’ predictions.

Relative to NAL, we see stronger dependence on radius for the disk mass models
from KF, DU, and RAL, which is particularly evident for the smaller chirp mass
case of Fig. 2. For the high chirp mass panel of Fig. 2, disagreement between
models is most pronounced in the highly compact domain of the parameter space.
The lower chirp mass panel demonstrates divergence between models in several
parts of the parameter space, notably in the top and bottom right corners. As
was the case for the dynamical ejecta, the differences between disk mass
models is at times comparable in magnitude to the ejecta predicted by the
models.

We then studied a series of 2396 best-fit equations of state derived from
aLIGO data for GW170817, and determined whether each equation of state could
reproduce some injected fiducial ejecta mass given the dynamical ejecta and
disk mass ejecta models we have considered. We find that there is slightly
stronger agreement between models for the Q=0.9 case, with some fiducial
observations consistently ruling out/accepting a subset of the proposed
equations of state. However, for most fiducial observations (and particularly
for more asymmetric systems), the large variations in the predictions of
different fitting formulae should urge us to proceed with extreme caution when
attempting to derive information about the equation of state of neutron stars
from kilonovae signals powered by BNS mergers.

For the case most similar to GW170817 ($M_{c}=1.186M_{\odot}$, Q=0.9,
$M_{tot}=0.05M_{\odot}$), we see all models agree that the least compact
(largest radius and, typically, largest maximum mass) equations of state are
consistent with such a system. This is one of the only spaces where all three
subfigures produce consistent results. Most models predict that the more
compact equations of state should be ruled out by such an observation (and
NAL, which does not at $Q=0.9$, would rule out compact neutron stars as well
if we had chosen a slightly more symmetric mass ratio).

Notably, we are unable to draw any conclusions about the NS EOS from the case
most similar to GW190425 ($M_{c}=1.44M_{\odot}$, $Q=0.72$, negligible ejecta).
As evidenced in Figure 5, each dynamical ejecta model provides a different
result. With KF’s model, the most compact EOSs are allowed and the lesser are
ruled out; DU’s model disallows all EOSs; and we see that no claims can be
made about the EOS when using NAL’s model. That the different models are able
to produce such discrepant results for unequal mass systems suggests that
further simulations in this regime are needed in order to properly calibrate
future ejecta models to a broader region of parameter space. We should also
note that this study is highly idealized, neglecting uncertainties due to the
ejecta composition and geometry, nuclear reactions, and photon transport in
the outflows. Our only objective is to assess (dis)agreement between existing
outflow models, and we see that even if differences between models were the
sole source of uncertainty, that uncertainty would already be extremely
limiting. However, with increasingly detailed simulations that cover a larger
set of parameters, we hope that in the future these models can be used to
place robust constraints on the amount of ejecta from NSNS mergers.

## Acknowledgements

A.H. and F.F. gratefully acknowledge support from the DOE through grant DE-
SC0020435, and from NASA through grant 80NSSC18K0565. GR and SMN are grateful
for financial support from the Nederlandse Organisatie voor Wetenschappelijk
Onderzoek (NWO) through the Projectruimte and VIDI grants (Nissanke). SMN also
acknowledges financial support from the NWO sector plan.

## References

  * Abbott  _et al._ [2017] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, and et al., Multi-messenger observations of a binary neutron star merger, The Astrophysical Journal 848, L12 (2017).
  * Wollaeger  _et al._ [2021] R. T. Wollaeger, C. L. Fryer, E. A. Chase, C. J. Fontes, M. Ristic, A. L. Hungerford, O. Korobkin, R. O’Shaughnessy, and A. M. Herring, A broad grid of 2d kilonova emission models (2021), arXiv:2105.11543 [astro-ph.HE] .
  * Barnes and Kasen [2013] J. Barnes and D. Kasen, EFFECT OF a HIGH OPACITY ON THE LIGHT CURVES OF RADIOACTIVELY POWERED TRANSIENTS FROM COMPACT OBJECT MERGERS, The Astrophysical Journal 775, 18 (2013).
  * Kasen  _et al._ [2017] D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. Ramirez-Ruiz, Origin of the heavy elements in binary neutron-star mergers from a gravitational-wave event, Nature 551, 80–84 (2017).
  * Bulla [2019] M. Bulla, possis: predicting spectra, light curves, and polarization for multidimensional models of supernovae and kilonovae, Monthly Notices of the Royal Astronomical Society 489, 5037–5045 (2019).
  * Shibata  _et al._ [2021] M. Shibata, S. Fujibayashi, and Y. Sekiguchi, Long-term evolution of neutron-star merger remnants in general relativistic resistive magnetohydrodynamics with a mean-field dynamo term, Phys. Rev. D 104, 063026 (2021), arXiv:2109.08732 [astro-ph.HE] .
  * Yamamoto  _et al._ [2008] T. Yamamoto, M. Shibata, and K. Taniguchi, Simulating coalescing compact binaries by a new code (sacra), Physical Review D 78, 10.1103/physrevd.78.064054 (2008).
  * Dietrich and Ujevic [2017] T. Dietrich and M. Ujevic, Modeling dynamical ejecta from binary neutron star mergers and implications for electromagnetic counterparts, Classical and Quantum Gravity 34, 105014 (2017).
  * Radice  _et al._ [2018] D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm, S. Bernuzzi, and L. F. Roberts, Binary neutron star mergers: Mass ejection, electromagnetic counterparts, and nucleosynthesis, The Astrophysical Journal 869, 130 (2018).
  * Krüger and Foucart [2020] C. J. Krüger and F. Foucart, Estimates for disk and ejecta masses produced in compact binary mergers, Physical Review D 101, 10.1103/physrevd.101.103002 (2020).
  * Nedora  _et al._ [2021] V. Nedora, F. Schianchi, S. Bernuzzi, D. Radice, B. Daszuta, A. Endrizzi, A. Perego, A. Prakash, and F. Zappa, Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers (2021), arXiv:2011.11110 [astro-ph.HE] .
  * Hotokezaka  _et al._ [2013] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Mass ejection from the merger of binary neutron stars, Physical Review D 87, 10.1103/physrevd.87.024001 (2013).
  * Bauswein  _et al._ [2013] A. Bauswein, S. Goriely, and H.-T. Janka, Systematics of dynamical mass ejection, nucleosynthesis, and radioactively powered electromagnetic signals from neutron-star mergers, The Astrophysical Journal 773, 78 (2013).
  * Dietrich  _et al._ [2015] T. Dietrich, S. Bernuzzi, M. Ujevic, and B. Brügmann, Numerical relativity simulations of neutron star merger remnants using conservative mesh refinement, Physical Review D 91, 10.1103/physrevd.91.124041 (2015).
  * Lehner  _et al._ [2016] L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero, E. O’Connor, M. Anderson, and D. Neilsen, Unequal mass binary neutron star mergers and multimessenger signals, Classical and Quantum Gravity 33, 184002 (2016).
  * Sekiguchi  _et al._ [2016] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, M. Shibata, and K. Taniguchi, Dynamical mass ejection from the merger of asymmetric binary neutron stars: Radiation-hydrodynamics study in general relativity, Physical Review D 93, 10.1103/physrevd.93.124046 (2016).
  * Dietrich  _et al._ [2017] T. Dietrich, M. Ujevic, W. Tichy, S. Bernuzzi, and B. Brügmann, Gravitational waves and mass ejecta from binary neutron star mergers: Effect of the mass ratio, Physical Review D 95, 10.1103/physrevd.95.024029 (2017).
  * Kiuchi  _et al._ [2019] K. Kiuchi, K. Kyutoku, M. Shibata, and K. Taniguchi, Revisiting the lower bound on tidal deformability derived by at 2017gfo, The Astrophysical Journal 876, L31 (2019).
  * Sekiguchi  _et al._ [2015] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Dynamical mass ejection from binary neutron star mergers: Radiation-hydrodynamics study in general relativity, Physical Review D 91, 10.1103/physrevd.91.064059 (2015).
  * Vincent  _et al._ [2020] T. Vincent, F. Foucart, M. D. Duez, R. Haas, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Unequal mass binary neutron star simulations with neutrino transport: Ejecta and neutrino emission, Physical Review D 101, 10.1103/physrevd.101.044053 (2020).
  * Perego  _et al._ [2019] A. Perego, S. Bernuzzi, and D. Radice, Thermodynamics conditions of matter in neutron star mergers, The European Physical Journal A 55, 10.1140/epja/i2019-12810-7 (2019).
  * Bernuzzi  _et al._ [2020] S. Bernuzzi, M. Breschi, B. Daszuta, A. Endrizzi, D. Logoteta, V. Nedora, A. Perego, D. Radice, F. Schianchi, F. Zappa, I. Bombaci, and N. Ortiz, Accretion-induced prompt black hole formation in asymmetric neutron star mergers, dynamical ejecta, and kilonova signals, Monthly Notices of the Royal Astronomical Society 497, 1488 (2020).
  * Nedora  _et al._ [2019] V. Nedora, S. Bernuzzi, D. Radice, A. Perego, A. Endrizzi, and N. Ortiz, Spiral-wave wind for the blue kilonova, The Astrophysical Journal 886, L30 (2019).
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimessenger constraints on the neutron-star equation of state and the hubble constant, Science 370, 1450–1453 (2020).
  * Hotokezaka  _et al._ [2011] K. Hotokezaka, K. Kyutoku, H. Okawa, M. Shibata, and K. Kiuchi, Binary neutron star mergers: Dependence on the nuclear equation of state, Phys. Rev. D 83, 124008 (2011).
  * Nedora  _et al._ [2020] V. Nedora, F. Schianchi, S. Bernuzzi, D. Radice, B. Daszuta, A. Endrizzi, A. Perego, A. Prakash, and F. Zappa, Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers (2020).
  * Coughlin  _et al._ [2019] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Multimessenger bayesian parameter inference of a binary neutron star merger, Monthly Notices of the Royal Astronomical Society: Letters 489, L91–L96 (2019).
  * Abbott  _et al._ [2020] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, and et al., Gw190425: Observation of a compact binary coalescence with total mass $\sim$ 3.4 $m_{\odot}$, The Astrophysical Journal 892, L3 (2020).
  * Camilletti  _et al._ [2022] A. Camilletti, L. Chiesa, G. Ricigliano, A. Perego, L. C. Lippold, S. Padamata, S. Bernuzzi, D. Radice, D. Logoteta, and F. M. Guercilena, Numerical relativity simulations of the neutron star merger GW190425: microphysics and mass ratio effects, arXiv (2022), arXiv:2204.05336 [astro-ph.HE] .
  * Chen  _et al._ [2021] H.-Y. Chen, S. Vitale, and F. Foucart, The Relative Contribution to Heavy Metals Production from Binary Neutron Star Mergers and Neutron Star–Black Hole Mergers, Astrophys. J. Lett. 920, L3 (2021), arXiv:2107.02714 [astro-ph.HE] .
  * Abbott  _et al._ [2018] B. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, V. Adya, and et al., Gw170817: Measurements of neutron star radii and equation of state, Physical Review Letters 121, 10.1103/physrevlett.121.161101 (2018).
  * Lindblom [2010] L. Lindblom, Spectral Representations of Neutron-Star Equations of State, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. Nissanke, F. Foucart, M. M. Kasliwal, M. Bulla, R. Fernandez, A. Henkel, T. Hinderer, K. Hotokezaka, K. Lukošiūtė, T. Venumadhav, S. Antier, M. W. Coughlin, T. Dietrich, and T. D. P. Edwards, The challenges ahead for multimessenger analyses of gravitational waves and kilonova: a case study on gw190425 (2021), arXiv:2102.11569 [astro-ph.HE] .


# Generation of bright collimated vortex $\gamma$-ray via laser driven cone-
fan target

Cui-Wen Zhang  Key Laboratory of Beam Technology of the Ministry of Education,
and College of Nuclear Science and Technology, Beijing Normal University,
Beijing 100875, China  Mamat-Ali Bake  School of Physics Science and
Technology, Xinjiang University, Urumqi 830046, China  Hong Xiao  Key
Laboratory of Beam Technology of the Ministry of Education, and College of
Nuclear Science and Technology, Beijing Normal University, Beijing 100875,
China  Hai-Bo Sang  Department of Physics, Beijing Normal University, Beijing
100875, China  Bai-Song Xie<EMAIL_ADDRESS>Key Laboratory of Beam
Technology of the Ministry of Education, and College of Nuclear Science and
Technology, Beijing Normal University, Beijing 100875, China  Institute of
Radiation Technology, Beijing Academy of Science and Technology, Beijing
100875, China

###### Abstract

We use numerical simulations to demonstrate that a source of bright collimated
vortex $\gamma$-ray with large orbital angular momentum can be achieved by
irradiating a circularly polarized laser with an intensity about
$10^{22}\rm{W/{cm^{2}}}$ on a cone-fan target. In the studied setup, electron
beam of energy of hundreds of MeV and vortex laser pulse are formed. And
furthermore a high quality vortex $\gamma$-ray is yielded with small
divergence of $5^{\circ}$ and high peak brilliance $\sim 5\times 10^{22}$
photons ${\rm\cdot s^{-1}\cdot mm^{-2}\cdot mrad^{-2}}$ $0.1\%\mathrm{BW}$ at
$10\mathrm{MeV}$. A considerable fraction of angular momentum of laser is
converted to electron beam and vortex $\gamma$-ray, which are roughly $27.8\%$
and $3\%$, respectively. And the conversion efficiency of energy from laser to
electron beam and vortex $\gamma$-ray are around $41\%$ and $3.8\%$. Moreover,
comparative simulations for different right radius of cone reveal that there
exists an optimal size that makes the highest angular momentum of $\gamma$-ray
photons to be around $2.8\times 10^{6}\hbar$. The comparative simulations for
different laser modes exhibit that it is more appropriate to choose the
circularly polarized laser to generate vortex $\gamma$-ray than the Laguerre-
Gaussian one.

###### pacs:

52.38.-r; 52.38.Ph; 52.65.Rr

##  I INTRODUCTION

The prospect of generating vortex $\gamma$-rays with orbital angular momentum
(OAM) has recently attracted particular interest due to its potential
applications in imaging Jesacher:2005 ; 2 , material science 3  and
astrophysics 4 ; 5 . Many vortex $\gamma$-rays generation schemes involving
circularly polarized (CP) lasers and/or Laguerre-Gaussian (LG) lasers have
been proposed 7 ; 8 ; 9 ; 10 ; 11 ; 12 ; 13 ; 14 ; 15 ; 16 ; 17 ;
channel2021HPL ; channel2022PST ; cone-foil2018APL ; 2018NJP . As we know, the
CP laser carries only spin angular momentum (SAM) of $n\delta\hbar$, where $n$
is the laser photon number, $\delta$ is the SAM quantum number of a laser
photon ($\delta=\pm 1$ for the right and left-hand) and $\hbar$ is the reduced
Planck constant. While the LG laser as a kind of vortex light, it has a
helical wave front, a twisted Poynting vector and hollow transverse field
distribution 6 , thus, it carries both of OAM and SAM as $n(l+\delta)\hbar$,
where $l$ represents the OAM quantum number of a laser photon, e.g., $l=0,\pm
1,\pm 2,\ldots$. Therefore, under the same intensity, the LG laser would carry
more angular momentum (AM), i.e., the sum of OAM and SAM, than the CP laser.
However, in a view point of practical experiment, the highest intensity of LG
laser achieved at present is around ${\rm 6.3\times 10^{19}W/cm^{2}}$
PRL2020LG , which is much less than that of the CP laser ${\rm
10^{23}W/cm^{2}}$ optica2021CP , thus, this strongly restricts the application
of the LG lasers for vortex $\gamma$-rays generation and related applications.

Vortex $\gamma$-rays can be generated by synchrotron radiation 7 ; 8  and
nonlinear Compton scattering (NCS) 9 ; 10 ; 11 ; 12 ; 13 ; 14 ; 15 ; 16 ; 17 ;
channel2021HPL ; channel2022PST ; cone-foil2018APL ; 2018NJP  processes in
quantum electrodynamics (QED) regime. The AM of laser can be transferred to
the OAM of electrons first, and then be further transferred to the OAM of
vortex $\gamma$-rays through synchrotron radiation. The AM conversion
efficiency is $1.43\%$ 7  and a vortex $\gamma$-ray with small divergence
angle of $10^{\circ}$ has been obtained 8 . In NCS, strong laser interacts
with counter-propagating energetic electron beam got by laser wakefield
acceleration 2018MRE  or laser direct acceleration 2020MRE , and vortex
$\gamma$-rays can be generated via $e^{-}+n\omega_{0}\rightarrow
e^{-}+\omega_{\gamma}$, where $n>1$ is an integer, $\omega_{0}$ and
$\omega_{\gamma}$ is the frequency of laser and vortex $\gamma$-ray,
respectively. In Ref. 10  Wang et al. got a vortex $\gamma$-ray with OAM of
$2.5\times 10^{18}\hbar$ and energy up to hundreds of $\rm{MeV}$. Chen et al.
16  pointed out that the AM of the absorbed laser photons is not solely
transferred to the emitted $\gamma$-ray photons, but due to radiation reaction
shared between the $\gamma$-ray photons and interacting electrons. In the
recent research Ref.channel2022PST . it is found that the AM conversion
efficiency is $2.1\%$. Therefore, NCS is a promising way for vortex
$\gamma$-rays generation based on currently available laser and accelerator
technologies. However, the enhancement of OAM, AM conversion efficiency,
energy and collimation degree of vortex $\gamma$-rays is still challenging and
need to be study furthermore.

To keep the advantage of CP laser intensity while overcome the deficiency of
AM, fortunately, previous studies revealed that light-fan, i.e., spiral-shaped
foil plasma PRL2014-fan ; channel2021HPL  is a useful plasma mirror
mirror2004PRE ; mirror2012NP ; mirror2013PRL  which can reflect the incident
CP laser and transfer it into a vortex one by adding phase $\exp(il\phi)$.
Compare to tube targets 11 ; channel2019PPCF ; channel2021HPL ; channel2022PST
, the cone targets can restrain laser better and the laser intensity can be
increased up to several tens of times, the laser oscillating field can pull
out the electrons in cone surface periodically and accelerate them forward via
laser pondermotive force and result in overdense energetic electron beams
cone2004POP ; cone2012POP ; cone2015POP ; cone2016JAP . Note that the AM
conversion efficiency from laser to vortex $\gamma$-ray is only $0.67\%$ in
the scheme by using tube target channel2021HPL . And the vortex $\gamma$-ray
in Ref.cone-foil2018APL  is found concentrating at $1\mathrm{MeV}$ when a CP-
LG laser pulse is irradiating on a cone-foil target. Thus, the cone-fan target
should be considered in vortex $\gamma$-ray generation.

Motivated by the cone-fan target advantages in both of increasing the laser
intensity, the produced electron energy and the high AM of reflected laser, in
this paper, we use three-dimensional (3D) particle-in-cell (PIC) simulations
to perform numerical study of vortex $\gamma$-ray generation via CP laser
irradiates a cone-fan target, note that both of the laser and fan are left-
hand in our study. The electrons on the inner surface of cone can be extracted
into cone-channel and then accelerated to the high energy. Moreover, the cone-
channel can guide these electrons effectively and result in overdense
energetic electron beam finally. Meanwhile, the SAM of CP laser is efficiently
converted into the OAM of electrons. When the incident laser reaches the fan
that attached on the rear side of cone, the laser is reflected and a vortex
one is formed. Then, the latter collides head-on with energetic electron beam
and a bright collimated vortex $\gamma$-ray is emitted through NCS. The OAM of
vortex $\gamma$-ray comes from both the vortex laser and electrons.

During the whole process, the fan is used as a reflecting mirror to get vortex
laser and spontaneously trigger the subsequent NCS. The cut-off energy and OAM
of vortex $\gamma$-ray are up to $200\mathrm{MeV}$ and $5.9\times
10^{18}\hbar$ respectively. Remarkably, the conversion rate of laser AM for
electron beam and vortex $\gamma$-ray are around $27.8\%$ and $3\%$, and
roughly $41\%$ of the laser energy is converted to electron beam and $3.8\%$
to vortex $\gamma$-ray. The peak brilliance of vortex $\gamma$-ray is about
$5\times 10^{22}$ photons ${\rm\cdot s^{-1}\cdot mm^{-2}\cdot mrad^{-2}}$
$0.1\%$ bandwidth at $10\mathrm{MeV}$ and divergence angle $5^{\circ}$. The
properties of this bright collimated vortex $\gamma$-ray might enable the
development of novel application in various domains.

Figure 1: (color online). The setup of the vortex $\gamma$-ray generation
scheme. A CP laser propagates forward along the $x$ axis from the left and
irradiates a cone-fan target. The electrons on the inner surface of the cone
can be extracted into the cone-channel and accelerated via direct laser
acceleration, and overdense energetic electron beam can be generated. The
incident laser is reflected when it reaches the fan and a vortex one is
formed. Then, the latter collides head-on with the energetic electron beam and
bright collimated vortex $\gamma$-ray is emitted through NCS.

##  II SCHEME SETUP AND SIMULATION PARAMETER

Figure 2: (color online). (a)-(d): Distributions of transverse electric field
$E_{y}$ of incident CP laser at different crossing sections at $t=15T_{0}$.
(e)-(h): Distributions of transverse electric field $E_{y}$ at different
crossing sections at $t=25T_{0}$ when the incident laser is completely
reflected by fan.

We performed 3D PIC simulations with open-source code EPOCH 2014epoch ;
2015epoch . As indicated in Fig.1, a left-hand Gaussian CP laser propagates
forward along the axis ($+x$ direction) of cone with a left-hand fan attached
behind. The dimensionless amplitude of CP laser field is
$a=a_{0}\exp(-(x/c-t_{0}))^{2}/\tau^{2})\exp({-(y^{2}+z^{2})/w^{2}})\left(\cos\varphi\widehat{\bm{e}}_{y}+\sin\varphi\widehat{\bm{e}}_{z}\right)$,
where $a_{0}=eE_{0}/m_{e}c\omega_{0}=100$ (corresponding to $I_{0}\approx
2.74\times 10^{22}\rm{W/{cm^{2}}}$) , $w=3\lambda_{0}$,
$\lambda_{0}=cT_{0}=1\mathrm{\mu m}$ and $\omega_{0}=2\pi/T_{0}$ are the
intensity, focal spot, wavelength and frequency of laser, respectively. $-e$
and $m_{e}$ are the electron charge and mass, $E_{0}$ is the electric field
amplitude, $c$ is the speed of light in vacuum, $\varphi$ is the laser phase
term and $\tau=10T_{0}$ is the laser pulse duration. The axial length of cone
is $15\lambda_{o}$ located from $x=1\lambda_{0}$ to $x=16\lambda_{0}$. The
inner left radius of cone $r_{1}=3\rm{\mu m}$, the right radius
$r_{2}=1.5\rm{\mu m}$ and the thickness is $1\lambda_{0}$. The fan consists of
eight parts with the same step height $\Delta=\lambda_{0}/16$ to mimic a
$\lambda_{0}/2$ spiral phase plate with thickness increasing with angle and
the minimum is $1\lambda_{0}$ to make sure that it can reflect the driving
laser, and the right end of the fan is located at $x=16\lambda_{0}$. The grid
size of simulation box is $40\lambda_{0}\times 12\lambda_{0}\times
12\lambda_{0}$ with $2000\times 300\times 300$ cells in the $x\times y\times
z$ direction and there are 9 macro-particles in each cell. The density of
electrons is $200n_{c}$, where $n_{c}=m_{e}\omega_{0}^{2}/4\pi e^{2}=1.1\times
10^{21}\rm{cm^{-3}}$ is the critical density.

In the first stage, shown in Fig.1, electrons within the inner surface of cone
are extracted into cone-channel periodically by oscillating laser field and
are accelerated forward via laser pondermotive force. Due to the guiding
action of cone, an overdense electron beam with energy of hundreds of
$\mathrm{MeV}$ is generated. During the process, the SAM of CP laser is
efficiently converted into the OAM of electrons, and we will discuss this
point later.

In the second stage, the incident CP laser reaches fan and is reflected into a
vortex one, see Fig.1. To illustrate this process, as shown in Fig.2, a
simulation without cone is performed. Figures 2(e)-2(h) indicate that the
reflected laser pulse is in vortex shape, and due to the lost in the
reflection process, the electric field of which is reduced to half that of the
incident laser. The laser mode of reflected pulse can be expanded by a series
of LG modes LGmode  and the amplitude of $\rm{LG_{nm}}$ can be defined by

$\displaystyle\rm{LG_{nm}}(r,\phi,x)=$
$\displaystyle(C_{nm}/w)\exp(-ikr^{2}/2R)\exp(-r^{2}/w^{2})$ (1)
$\displaystyle\times\exp[-i(n+m+1)\psi]$
$\displaystyle\times\exp[-i(n-m)\phi](-1)^{p}$
$\displaystyle\times(r\sqrt{2}/w)^{|l|}L_{p}^{|l|}(2r^{2}/w^{2}),$

where

$\displaystyle R(x)=(x_{R}^{2}+x^{2})/x,$ (2)

$\displaystyle w(x)=[2(x_{R}^{2}+x^{2})/kx_{R}]^{1/2},$ (3)

$\displaystyle\psi(x)=\arctan(x/x_{R}),$ (4)

$r$, $\phi$ and $x$ are cylindrical coordinates, $C_{nm}$ is the normalization
constant, $k=2\pi/\lambda$ is the wave number, $x_{R}$ is the Rayleigh range,
${L_{p}}^{l}(2r^{2}/w^{2})$ is the generalized Laguerre polynomial, $l=|n-m|$
is the azimuthal mode index, and $p=\min(n,m)$ is the number of radial nodes.
The Gaussian mode of incident laser beam in simulation is $\rm{LG_{00}}$ whose
wave front has been modified by fan structure and the mode decomposition can
be written using expansion coefficients LGmode

$\displaystyle a_{st}=$
$\displaystyle\left\langle\mathrm{LG}_{st}|\exp(-i\Delta\phi)|\mathrm{LG}_{nm}\right\rangle$
(5) $\displaystyle=$ $\displaystyle\iint
rdrd\phi\left(C_{st}^{*}/w_{st}\right)\exp\left(ikr^{2}/2R_{st}-r^{2}/w_{st}^{2}\right)$
$\displaystyle\times\exp[-i(s-t)\phi](-1)^{\min(s,t)}\left(r\sqrt{2}/w_{st}\right)^{|s-t|}$
$\displaystyle\times
L_{\min(s,t)}^{|s-t|}\left(2r^{2}/w_{st}^{2}\right)\exp(-i\Delta\phi)$
$\displaystyle\times\left(C_{00}/w\right)\exp\left(-ikr^{2}/2R-r^{2}/w^{2}\right),$

where

$\Delta\phi=\sum_{n=0}^{7}H\left(\phi-\frac{\pi}{4}n\right)H\left(\frac{\pi}{4}n+\frac{\pi}{4}-\phi\right)\frac{2\pi}{8}n,$
(6)

and $H(x)$ is the Heaviside function. As $\Delta\phi\approx\phi$, only modes
with $(s-t)-(n-m)=1$ contribute most in the $\phi$ integral and the percent of
modes is given by $I_{st}=|a_{st}|^{2}$. Rayleigh range and waist are assumed
to be equal to those of incident laser. Our calculations show that
$I_{10}\approx 79.4\%$, $I_{21}\approx 9.9\%$ and $I_{32}\approx 3.8\%$, thus,
the dominant mode of reflected laser pulse is $\rm{LG_{10}}$ which is
consistent with the simulation result in Fig.2.

In the third stage, reflected vortex laser collides head-on with energetic
electron beam and copious hundreds of $\mathrm{MeV}$ $\gamma$-ray photons are
emitted through NCS, see Fig.1. The QED parameter
$\chi_{e}=e\hbar(m^{3}_{e}c^{4})|F_{\mu\nu}p^{\nu}|=(\gamma_{e}/E_{s})\sqrt{(\textbf{E}+\bm{\beta}\times\textbf{B})^{2}-(\bm{\beta}\cdot\textbf{E})^{2}}$
dominants photon emission probability QEDparameter1985 , where $F_{\mu\nu}$ is
the electromagnetic field tensor, $p^{\nu}$ is the electron four-momentum and
$E_{s}=(m_{e}^{2}c^{3})/e\hbar\approx 1.3\times 10^{18}\rm{Vm^{-1}}$
(corresponding to $I_{s}\approx 4.6\times 10^{29}\rm{W/{cm^{2}}}$) is the
Schwinger electric field, at that intensity electrons can be accelerated to
$m_{e}c^{2}$ in Compton wavelength distance
$\lambda_{c}=\hbar/m_{e}c=3.8616\times 10^{-13}\rm{m}$ and the vacuum decays
to create pairs Schwinger1951 . $\bm{\beta}=\textbf{v}/c$ is the normalized
velocity of electron, $\gamma_{e}$ is the electron relativistic factor, and E
and B are the electric and magnetic fields acting on electrons. It is worth
noting that radiation barely occurs when energetic electron beam co-propagates
with intense laser in which $\chi_{e}\rightarrow 0$. While, when electron beam
and laser counter-propagate, $\chi_{e}$ can be described as
$\chi_{e}=2\gamma_{e}|E_{\perp}|/E_{s}$ and the value of which is large and
copious $\gamma$-ray photons can be emitted. During the process of NCS, the AM
of vortex laser and electron beam are converted to $\gamma$-ray photons, and
we will discuss this point later.

##  III RESULTS AND DISCUSSION

###  III.1 Generation of vortex electron beam and $\gamma$-ray

Figure 3: (color online). (a) Three-dimensional isosurface distribution of
electron beam at $t=15T_{0}$. (b) Energy spectrum of electrons at $t=10T_{0}$,
$20T_{0}$, $25T_{0}$ and $30T_{0}$.

As indicated in Fig.3(a), energetic electron beam appears in spiral shape,
which results from the helical structure of laser electric field. Figure 3(b)
presents the energy spectrum of the electrons at different time. Under the
action of the ponderomotive force of laser, the cut-off energy of electrons
increases and reaches a maximum of $\sim\rm{500MeV}$ at $20T_{0}$. And after
which, it starts to decrease due to the emission of vortex $\gamma$-ray. The
main emission process stopped at about $25T_{0}$ when CP laser pulse is
completely reflected, and the cut-off energy of electrons no longer changes
significantly.

Figure 4: (color online). (a) Energy spectrum of $\gamma$-ray photons at
$t=18T_{0}$, $20T_{0}$, $22T_{0}$ and $30T_{0}$. (b) Three-dimensional
isosurface distribution of vortex $\gamma$-ray at $t=20T_{0}$.

$\gamma$-ray photons are emitted via NCS when the reflected vortex laser and
energetic electron beam collides head-on. As shown in Fig.4(a), the cut-off
energy of $\gamma$-ray photons increases and reaches a maximum of
$\sim\rm{200MeV}$ at about $22T_{0}$. It is worth noting that the isosurface
of vortex $\gamma$-ray, shown in Fig.4(b), is also helical. Which is because
the direction of instantaneously emitted $\gamma$-ray photon is at an angle
$1/\gamma_{e}$ with that of moving electron, and it is negligible for
ultrarelativistic electrons ($\gamma_{e}>>1$). Besides, the $\chi_{e}$ at
$20T_{0}$ is $0.93$, which coincides with the calculation result by using
$\hbar\omega_{\gamma}\approx 0.44\chi_{e}\gamma m_{e}c^{2}$.

Figure 5: (color online). Evolution of $M_{l}$, $M_{e^{-}}$, $M_{\gamma}$ and
$M_{\gamma-ph}$ for the AM of laser (black), electron beam (red), vortex
$\gamma$-ray (blue) and average of $\gamma$-ray photons (green).

In order to figure out the evolution of AM over time, we tracked the different
AM, i.e., $M_{l}$, $M_{e^{-}}$, $M_{\gamma}$ and $M_{\gamma-ph}$ for the
laser, electron beam, vortex $\gamma$-ray and average of $\gamma$-ray photons,
as shown in Fig.5. The AM of Gaussian CP laser can be calculated by

${M}_{l}=\big{|}\varepsilon_{0}\int\textbf{r}\times(\textbf{E}\times\textbf{B})\rm{d\nu}\big{|}=\big{|}\varepsilon_{0}\triangle\nu\sum_{i}\textbf{r}_{i}\times(\textbf{E}_{i}\times\textbf{B}_{i})\big{|},$
(7)

where $\varepsilon_{0}$ is the permittivity of vacuum, r is the grid position,
$\triangle\nu$ is the volume of one cell and $i$ is the ordinal number of
cell. So, the averaged AM of each laser photon is

$M_{l-ph}=\frac{\big{|}\varepsilon_{0}\int\textbf{r}\times(\textbf{E}\times\textbf{B})\rm{d\nu}\big{|}}{\frac{1}{2}\int(\varepsilon_{0}\textbf{E}^{2}+\frac{1}{\mu_{0}}\textbf{B}^{2})\rm{d\nu}}\hbar\omega_{0}=(l+\delta)\hbar,$
(8)

where
$\frac{1}{2}\int(\varepsilon_{0}\textbf{E}^{2}+\frac{1}{\mu_{0}}\textbf{B}^{2})\rm{d\nu}$
is the total electromagnetic energy of laser, $\mu_{0}$ is the vacuum
permeability. The black line in Fig.5 shows that the $M_{l}$ reaches the
maximum of $1.58\times 10^{20}\hbar$ which is close to that in the vacuum
$2\times 10^{20}\hbar$ at $10T_{0}$ when the CP laser completely enters the
simulation region, and the difference is transferred to the electrons in cone.

The AM of particle beam with respect to the $x$-axis is

$M_{particle}=\big{|}\sum_{j}(\textbf{r}_{j}\times\textbf{p}_{j})\big{|}=\big{|}\sum_{j}(yp_{z}-zp_{y})_{j}\hat{\textbf{x}}\big{|},$
(9)

where $\textbf{r}_{j}$ and $\textbf{p}_{j}$ are the position and momentum of
the $j$-th particle respectively, and $p_{y}$ and $p_{z}$ are the momentum in
the $y$ and $z$ direction. As shown by the red line in Fig.5, the $M_{e^{-}}$
increases since CP laser propagates into the cone-channel and reaches maximum
of $5.56\times 10^{19}\hbar$ at $23T_{0}$ and then decreases. In the first
stage, before $15T_{0}$, electrons are accelerated by CP laser and the SAM of
laser is converted to the OAM of electrons. In the second stage, the electrons
are still being accelerated because the intensity of the forward laser is
greater than that of the reflected laser, thus, $M_{e^{-}}$ still increases.
In the third stage, due to radiation reaction, the AM of the reflected laser
can both be transferred to the electrons and $\gamma$-ray photons 16 .
Therefore, the $M_{e^{-}}$ continues to increase. After the third stage, i.e.,
$23T_{0}$, there are no sources of OAM for electrons, and $M_{e^{-}}$
decreases gradually.

Figure 6: (color online). (a) Evolution of conversion efficiency of laser AM
for electron beam $\eta^{M}_{e^{-}}$ (black) and vortex $\gamma$-ray
$\eta^{M}_{\gamma}$ (red). (b) Evolution of conversion efficiency of laser
energy for electron beam $\eta^{\varepsilon}_{e^{-}}$ (black) and vortex
$\gamma$-ray $\eta^{\varepsilon}_{\gamma}$ (red).

As shown by the blue line of $M_{\gamma}$ and the green line of $M_{\gamma-
ph}$ in Fig.5, both of them increase since $t=15T_{0}$ when NCS occurs and
reach maximums of $5.9\times 10^{18}\hbar$ at $26T_{0}$ and $2.8\times
10^{6}\hbar$ at $20T_{0}$, respectively. And the OAM of $\gamma$-ray photons
comes from that of the reflected laser. It is worth noting that $M_{\gamma}$
and $M_{\gamma-ph}$ reach maximums at different time, this is because there
are always $\gamma$-ray photons being emitted after $15T_{0}$ and the number
of $\gamma$-ray photons at $26T_{0}$ is more than that at $20T_{0}$.

Figure 7: (color online). (a) The polar-energy distribution of $\gamma$-ray
photons at $t=20T_{0}$. The red line exhibits polar-number distribution. (b)
The angular distribution of $\gamma$-ray photons at $t=20T_{0}$. (c) The
position distribution of $\gamma$-ray photons at $t=20T_{0}$.

The conversion efficiency of laser AM for particles can be calculated by

$\eta^{M}_{particle}=\frac{|\sum_{j}(\textbf{r}_{j}\times\textbf{p}_{j})|}{|\varepsilon_{0}\triangle\nu\sum_{i}\textbf{r}_{i}\times(\textbf{E}_{i}\times\textbf{B}_{i})|}.$
(10)

So, the conversion efficiencies of laser AM to electron beam
$\eta^{M}_{e^{-}}$ and vortex $\gamma$-ray $\eta^{M}_{\gamma}$ are roughly
$27.8\%$ and $3\%$, as shown in Fig.6(a), and $\eta^{M}_{\gamma}$ is almost 4
times higher than that of the channel scheme in Ref.channel2021HPL  and 1.5
times higher than that of Ref.channel2022PST . Furthermore, at the end of the
simulation, the laser energy conversion efficiencies to electron beam
$\eta^{\varepsilon}_{e^{-}}$ and vortex $\gamma$-ray
$\eta^{\varepsilon}_{\gamma}$ are around $41\%$ and $3.8\%$, see Fig.6(b),
which are double those in Ref.channel2021HPL .

The polar angle of particle can be calculated by
$\arcsin(\sqrt{p^{2}_{\bot}/p^{2}})$ for $p_{x}\geq 0$ and
$\pi-\arcsin(\sqrt{p^{2}_{\bot}/p^{2}})$ for $p_{x}<0$, where
$p^{2}_{\bot}=p^{2}_{y}+p^{2}_{z}$ and $p^{2}=p^{2}_{\bot}+p^{2}_{x}$. Figure
7(a) indicates that $\gamma$-ray photons are concentrated in a narrow range of
$5^{\circ}$ centered around polar angle $7^{\circ}$, the corresponding photon
energy is $10\rm{MeV}$ and the number of $\gamma$-ray photons is about
$1\times 10^{9}$, which are all rather better. Figures 7(b) and 7(c) show that
the distribution of momentum and position of $\gamma$-ray photons in $y-z$
plane are in the shape of a doughnut, which is consistent with the helical
structure of the vortex $\gamma$-ray. From Fig.7(c), it can be roughly
estimated that $\gamma$-ray photons are distributed in the radial range
$1\lambda_{0}$-$1.2\lambda_{0}$. And as can be seen from Fig.4(b), the
duration of vortex $\gamma$-ray pulse $\tau_{1}\approx
6\lambda_{0}/c=20\rm{fs}$. Thus, the peak brilliance of vortex $\gamma$-ray is
$\sim 5\times 10^{22}$ photons ${\rm\cdot s^{-1}\cdot mm^{-2}\cdot mrad^{-2}}$
$0.1\%$ bandwidth at $10\rm{MeV}$, which is promising for some potential
applications in various domains.

###  III.2 Vortex $\gamma$-ray emission of different schemes

Figure 8: (color online). The evolution of $M_{\gamma}$ (red), $\eta^{\varepsilon}_{\gamma}$ (blue) and $M_{\gamma-ph}$ (black) with the increase of right radius $r_{2}$.  Table 1: The maximums of $M_{l}$, $M_{e^{-}}$, $M_{\gamma}$, $M_{\gamma-ph}$, $\eta^{M}_{e^{-}}$ and $\eta^{M}_{\gamma}$ of three cases of CP-cone-fan, LG-cone-fan and LG-cone-foil. | case | $M_{l}(\times 10^{18}\hbar)$ | $M_{e^{-}}(\times 10^{18}\hbar)$ | $M_{\gamma}(\times 10^{18}\hbar)$ | $M_{\gamma-ph}(\times 10^{6}\hbar)$ | $\eta^{M}_{e^{-}}$ | $\eta^{M}_{\gamma}$
---|---|---|---|---|---|---|---
| CP-cone-fan | $-200$ | $-55.6$ | $-5.9$ | $-2.8$ | $27.8$ | $2.9$
| LG-cone-fan | $-400$ | $-92.9$ | $-2.0$ | $-1.1$ | $23.2$ | $0.5$
| LG-cone-foil | $-400$ | $-98.6$ | $-2.1$ | $-0.9$ | $24.7$ | $0.5$
Figure 9: (color online). (a)-(d): Distributions of transverse electric field
$E_{y}$ of incident $\rm{LG_{10}}$ laser at different crossing sections at
$t=15T_{0}$. (e)-(h): Distributions of transverse electric field $E_{y}$ at
different crossing sections at $t=25T_{0}$ when the incident laser is
completely reflected by fan.

We have performed comparative simulations for several different right radius
$r_{2}$ while keeping other parameters the same. Figure 8 indicates that the
red line of $M_{\gamma}$ and the blue line of $\eta^{\varepsilon}_{\gamma}$
all decrease with the increase of $r_{2}$. As $r_{2}$ increases, the cone
becomes less constrained to the incident laser and fewer electrons are dragged
out of the cone. Therefore, the average energy and density of the electrons
decrease cone2016JAP . And accordingly, the $M_{\gamma}$ and
$\eta^{\varepsilon}_{\gamma}$ decrease. However, $M_{\gamma-ph}$, the black
line in Fig.8, reaches maximum of $\sim 2.8\times 10^{6}\hbar$ at
$r_{2}=1.5\rm{\mu m}$ which indicates that $r_{2}=1.5\rm{\mu m}$ is the best
choice when we want the highest AM of $\gamma$-ray photons.

We also studied the case of CP-LG laser, which remained the same laser energy
as the CP one. The LG laser is left-hand and in $\rm{LG_{10}}$ mode, and the
target attached behind are fan and foil respectively, i.e., LG-cone-fan and
LG-cone-foil cases. The equation of $\rm{LG_{10}}$ laser is

$\displaystyle\mathrm{LG}_{10}=$
$\displaystyle(C_{10}/w)\exp(-ikr^{2}/2R)\exp\left({-r^{2}}/{w^{2}}\right)$
(11) $\displaystyle\times\exp(-i\psi)\exp(-i\phi)$
$\displaystyle\times\left({\sqrt{2}r}/{w}\right)^{1}L_{0}^{1}\left({2r^{2}}/{w^{2}}\right).$

So, equation 5 becomes

$\displaystyle a_{st}=$
$\displaystyle\left\langle\mathrm{LG}_{st}|\exp(-i\Delta\phi)|\mathrm{LG}_{nm}\right\rangle$
(12) $\displaystyle=$ $\displaystyle\iint
rdrd\phi\left(C_{st}^{*}/w_{st}\right)\exp\left(ikr^{2}/2R_{st}-r^{2}/w_{st}^{2}\right)$
$\displaystyle\times\exp[-i(s-t)\phi](-1)^{\min(s,t)}\left(r\sqrt{2}/w_{st}\right)^{|s-t|}$
$\displaystyle\times
L_{\min(s,t)}^{|s-t|}\left(2r^{2}/w_{st}^{2}\right)\exp(-i\Delta\phi)$
$\displaystyle\times\left(C_{10}/w\right)\exp\left(-ikr^{2}/2R-r^{2}/w^{2}\right)$
$\displaystyle\times\exp(-i\phi)\left(r\sqrt{2}/w\right)^{1}\times
L_{0}^{1}\left(2r^{2}/w^{2}\right),$

where $L_{0}^{1}\left(\frac{2r^{2}}{w^{2}}\right)=1$. The photon of
$\rm{LG_{10}}$ laser carries both OAM and SAM of $(l+\delta)\hbar=-2\hbar$,
which is twice that of CP laser, see Table 1. Consequently, the $M_{e^{-}}$
and $\eta^{M}_{e^{-}}$ in both cases are almost twice that in CP case. During
reflection, modes with $(s-t)-(n-m)=1$ contribute most, and our calculations
show that $I_{20}\approx 89\%$ and $I_{31}\approx 7.4\%$. While, $M_{\gamma}$,
$M_{\gamma-ph}$ and $\eta^{M}_{\gamma}$ are all lower than those in CP case
which is because that the transverse electric field $E_{y}$ of the reflected
laser is much lower than that of incident $\rm{LG_{10}}$ laser, as shown in
Fig.9. And the different is lost during the reflection process. Therefore, the
CP laser is much more appropriate than LG laser in our scheme. Additionally,
comparative simulations also indicate that the AM of vortex $\gamma$-ray is
influenced by the mode of laser and the number of steps of fan.

##  IV CONCLUSION

In conclusion, we have proposed an all-optical scheme to generate a bright
collimated vortex $\gamma$-ray with large OAM by simulating a CP laser
irradiates on a cone-fan target. In contrast to other schemes for generating
vortex $\gamma$-rays, the SAM of CP laser is converted into the OAM of
electron beam and vortex $\gamma$-ray with high efficiencies of $27.8\%$ and
$3\%$, the energy conversion efficiency of laser for electron beam and vortex
$\gamma$-ray are around $41\%$ and $3.8\%$, and the OAM of vortex $\gamma$-ray
and averaged $\gamma$-ray photons are $6\times 10^{18}\hbar$ and $3\times
10^{6}\hbar$. The peak brilliance of vortex $\gamma$-ray is about $5\times
10^{22}$ photons ${\rm\cdot s^{-1}\cdot mm^{-2}\cdot mrad^{-2}}$ $0.1\%$
bandwidth at $10\rm{MeV}$ and divergence angle $5^{\circ}$. The properties of
this bright collimated vortex $\gamma$-ray might enable the development of
novel application in materials and laboratory astrophysics.

Comparative simulations of different right radius $r_{2}$ show that
$M_{\gamma-ph}$ reaches maximum at $r_{2}=1.5\rm{\mu m}$, and comparative
simulations of different laser modes indicate that CP laser is better than LG
laser in our scheme. In addition, it can be inferred from our study that both
the mode of laser and the number of steps of fan can also affect the
generation of vortex $\gamma$-rays, which is worthy of further study.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(NSFC) under Grant No.11875007 and No. 11935008. The computation was carried
out at the HSCC of the Beijing Normal University. We acknowledgment the open
source PIC code EPOCH.

## References

  * (1) S. F$\mathrm{\ddot{u}}$rhapter, A. Jesacher, S. Bernet, and M. Ritsch-Marte, Spiral interferometry, Opt. Lett. 30, 1953 (2005).
  * (2) C. Maurer, A. Jesacher, S. Bernet, and M. Ritsch-Marte, What spatial light modulators can do for optical microscopy, Laser Photon. Rev. 5, 81-101 (2010).
  * (3) K. W. D. Ledingham, and W. Galster, Laser-driven particle and photon beams and some applications, New J. Phys. 12, 045005 (2010).
  * (4) T. Takiwaki, K. Kotake, and K. Sato, Special relativistic simulations of magnetically dominated jets in collapsing massive stars, Astrophys. J. 691, 1360 (2009).
  * (5) F. Tamburini, B. Thid$\mathrm{\acute{e}}$, G. Molina-Terriza, and G. Anzolin, Twisting of light around rotating black holes, Nat. Phys. 7, 195 (2011).
  * (6) C. Liu, B. F. Shen, X. M. Zhang, Y. Shi, L. L. Ji, W. P. Wang, L. Q. Yi, L. G. Zhang, T. J. Xu, Z. K. Pei, and Z. Z. Xu, Generation of $\gamma$-ray beam with orbital angular momentum in the QED regime, Phys. Plasmas 23, 093120 (2016).
  * (7) L. B. Ju, C. T. Zhou, T. W. Huang, K. Jiang, C. N. Wu, T. Y. Long, L. Li, H. Zhang, M. Y. Yu, and S. C. Ruan, Generation of Collimated Bright Gamma Rays with Controllable Angular Momentum Using Intense Laguerre-Gaussian Laser Pulses, Phys. Rev. A 12, 014054 (2019)
  * (8) S. Stock, A. Surzhykov, S. Fritzsche, and D. Seipt, Compton scattering of twisted light: Angular distribution and polarization of scattered photons, Phys. Rev. A 92, 013401 (2015).
  * (9) Z. Gong , R. H. Hu, H. Y. Lu, J. Q. Yu, D. H. Wang, E. G. Fu, C. E. Chen, X. T. He, and X. Q. Yan, Brilliant GeV gamma-ray flash from inverse Compton scattering in the QED regime, Plasma Phys. Control. Fusion 60, 044004 (2018).
  * (10) B. Feng, C. Y. Qin, X. S. Geng, Q. Yu, W. Q. Wang, Y. T. Wu, X. Yan, L. L. Ji, and B. F. Shen, The emission of $\gamma$-Ray beams with orbital angular momentum in laser-driven micro-channel plasma target, Sci. Rep. 9, 18780 (2019).
  * (11) J. Wang, X. B. Li, L. F. Gan, Y. Xie, C. L. Zhong, C. T. Zhou, S. P. Zhu, X. T. He, and B. Qiao, Generation of Intense Vortex Gamma Rays via Spin-to-Orbital Conversion of Angular Momentum in Relativistic Laser-Plasma Interactions, Phys. Rev. A 14, 014094 (2020).
  * (12) Y. Taira, T. Hayakawa, and M. Katoh, Gamma-ray vortices from nonlinear inverse Thomson scattering of circularly polarized light, Sci. Rep. 7, 5018 (2017).
  * (13) Y. Taira, and M. Katoh, Gamma-ray vortices emitted from nonlinear inverse Thomson scattering of a two-wavelength laser beam, Phys. Rev. A 98, 052130 (2018).
  * (14) V. Petrillo, G. Dattoli, I. Drebot, and F. Nguyen, Compton Scattered X-Gamma Rays with Orbital Momentum, Phys. Rev. Lett. 117, 123903 (2016).
  * (15) Y. Y. Chen, J. X. Li, K. Z. Hatsagortsyan, and C. H. Keitel, $\gamma$-Ray Beams with Large Orbital Angular Momentum via Nonlinear Compton Scattering with Radiation Reaction, Phys. Rev. Lett. 121, 074801 (2018).
  * (16) Y. Y. Chen, K. Z. Hatsagortsyan, and C. H. Keitel, Generation of twisted $\gamma$-ray radiation by nonlinear Thomson scattering of twisted light, Matter Radiat. Extremes 4, 024401 (2019).
  * (17) H. Zhang, J. Zhao, Y. T. Hu, Q. N. Li, Y. Lu, Y. Cao, D. B. Zou, Z. M. Sheng, F. Pegoraro, P. McKenna, F. Q. Shao, and T. P. Yu, Efficient bright $\gamma$-ray vortex emission from a laser-illuminated light-fan-in-channel target, High Power Laser Sci. Eng. 9, 1-11 (2021).
  * (18) M. A. Bake, S. Tang, and B. S. Xie, Bright $\gamma$-ray source with large orbital angular momentum from the laser near-critical-plasma interaction, Plasma Sci. Technol. 24, 095001 (2022).
  * (19) X. L. Zhu, M. Chen, T. Pu. Yu, S. M. Weng, L. X. Hu, P. McKenna, and Z. M. Sheng, Bright attosecond $\gamma$-ray pulses from nonlinear Compton scattering with laser illuminated compound targets, Appl. Phys. Lett. 112, 174102 (2018).
  * (20) X. L. Zhu, T. P. Yu, M. Chen, S. M. Weng, and Z. M. Sheng, Generation of GeV positron and $\gamma$-photon beams with controllable angular momentum by intense lasers, New J. Phys. 20, 083013 (2018).
  * (21) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes, Phys. Rev. A 45, 8185 (1992).
  * (22) W. P. Wang, C. Jiang, H. Dong, X. M. Lu, J. F. Li, R. J. Xu, Y. J. Sun, L. H. Yu, Z. Guo, X. Y. Liang, Y. X. Leng, R. X. Li, and Z. Z. Xu, Hollow Plasma Acceleration Driven by a Relativistic Reflected Hollow Laser, Phys. Rev. Lett. 125, 034801 (2020).
  * (23) J. W. Yoon, Y. G. Kim, I. W. Choi, J. H. Sung, H. W. Lee, S. K. Lee, and C. H. Nam, Realization of laser intensity over $10^{23}\rm{W/cm^{2}}$, Optica 8, 630-635 (2021).
  * (24) Y. Shi, B. F. Shen, L. G. Zhang, X. M. Zhang, W. P. Wang, and Z. Z. Xu, Light Fan Driven by a Relativistic Laser Pulse, Phys. Rev. Lett. 112, 235001 (2014).
  * (25) G. Doumy, F. Qu$\mathrm{\acute{e}}$r$\mathrm{\acute{e}}$, O. Gobert, M. Perdrix, Ph. Martin, P. Audebert, J. C. Gauthier, J. P. Geindre, and T. Wittmann, Complete characterization of a plasma mirror for the production of high-contrast ultraintense laser pulses, Phys. Rev. E 69, 026402 (2004).
  * (26) K. T. Phuoc, S. Corde, C. Thaury, V. Malka, A. Tafzi, J. P. Goddet, R. C. Shah, S. Sebban, and A. Rousse, All-optical Compton gamma-ray source, Nat. Photon 6, 308-311 (2012).
  * (27) S. Chen, N. D. Powers, I. Ghebregziabher, C. M. Maharjan, C. Liu, G. Golovin, S. Banerjee, J. Zhang, N. Cunningham, A. Moorti, S. Clarke, S. Pozzi, and D. P. Umstadter, $\rm{MeV}$-Energy X Rays from Inverse Compton Scattering with Laser-Wakefield Accelerated Electrons, Phys. Rev. Lett. 110, 155003 (2013).
  * (28) M. H. Li, L. M. Chen, D. Z. Li, K. Huang, Y. F. Li, Y. Ma, W. C. Yan, M. Z. Tao, J. H. Tan, Z. M. Sheng, and J. Zhang, Collimated gamma rays from laser wakefield accelerated electrons, Matter Radiat. Extremes 3, 188-196 (2018).
  * (29) K. Xue, Z. K. Dou, F. Wan, T. P. Yu, W. M. Wang, J. R. Ren, Q. Zhao, Y. T. Zhao, Z. F. Xu, and J. X. Li, Generation of highly-polarized high-energy brilliant $\gamma$-rays via laser-plasma interaction, Matter Radiat. Extremes 5, 054402 (2020).
  * (30) T. Y. Long, C. T. Zhou, T. W. Huang, K. Jiang, L. B. Ju, H. Zhang, T. X. Cai, M. Y. Yu, B. Qiao, S. C. Ruan, and X. T. He, All-optical generation of petawatt gamma radiation via inverse Compton scattering from laser interaction with tube target, Plasma Phys. Control. Fusion 61, 085002 (2019).
  * (31) Y. Sentoku, K. Mima, H. Ruhl, Y. Toyama, R. Kodama, and T. E. Cowan, Laser light and hot electron micro focusing using a conical target, Phys. Plasmas 11, 3083 (2004).
  * (32) F. Liu , X. X. Liu , B. C. Liu , W. J. Ding, F. Du, Y. T. Li, J. L. Ma, and X. L. Liu et al., Micro focusing of fast electrons with opened cone targets, Phys. Plasmas 19, 045004 (2012).
  * (33) L. X. Hu, T. P. Yu, F. Q. Shao, D. B. Zou, and Y. Yin, Enhanced dense attosecond electron bunch generation by irradiating an intense laser on a cone target, Phys. Plasmas 22, 033104 (2015).
  * (34) L. X. Hu, T. P. Yu, F. Q. Shao, W. Luo, and Y. Yin, A bright attosecond X-ray pulse train generation in a double-laser-driven cone target, J. Appl. Phys. 119, 243301 (2016).
  * (35) C. P. Ridgers, J. G. Kirk, R. Duclous, T. G. Blackburn, C. S. Brady, K. Bennett, T. D. Arber, and A. R. Bell, Modelling gamma-ray photon emission and pair production in high-intensity laser-matter interactions, J. Comput. Phys. 260, 273-285 (2014).
  * (36) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, Contemporary particle-in-cell approach to laser-plasma modelling, Plasma Phys. Control. Fusion 57, 113001 (2015).
  * (37) M. W. Beijersbergen, R. P. C. Coerwinkel, M. Kristensen, and J. P. Woerdman, Helical-wavefront laser beams produced with a spiral phaseplate, Opt. Commun. 112(5-6), 321-327 (1994).
  * (38) V. I. Ritus, Quantum Effects of the Interaction of Elementary Particles with an Intense Electromagnetic Field, J. Sov. Laser Res. 6(5), 497-617 (1985).
  * (39) J. Schwinger, On Gauge Invariance and Vacuum Polarization, Phys. Rev. 82, 664-679 (1951).