# Impact of Asymmetrical Mass Ejection from Proto-White Dwarfs on the
Properties of Binary Millisecond Pulsars

Wen-Shi Tang${}^{1,~{}2}$, Shi-Jie Gao${}^{1,~{}2}$, Xiang-Dong
Li${}^{1,~{}2}$
1School of Astronomy and Space Science, Nanjing University, Nanjing, 210023,
People’s Republic of China
2Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University,
Ministry of Education, Nanjing, 210023, People’s Republic of China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The standard formation theory of binary millisecond pulsars (BMSPs) predicts
efficient orbital circularization due to tidal interaction during the previous
mass transfer phase. Therefore, BMSPs are expected to have a circular orbit.
However, the discovery of several eccentric BMSPs (eBMSPs) with a white dwarf
(WD) companion has challenged this picture. In particular, recent observation
reveals that the spin angular momentum of the eBMSP J0955$-$6150 is tilted at
an angle $>4.8^{\rm\degree}$ from the orbital angular momentum. This is the
first time that a tilt angle is deduced for eBMSPs, which provides an
important clue to their formation mechanism. Both the orbital eccentricity and
tilt angle could be qualitatively accounted for by asymmetrical mass ejection
during thermonuclear flashes from proto-WDs (so-called the thermonuclear
rocket model), but detailed studies are still lacking. In this paper, we
simulate the impact of the kick caused by asymmetrical mass ejection on the
properties of BMSPs. We find that the thermonuclear rocket model can
potentially explain the observational characteristics of both eBMSPs and
normal BMSPs under reasonable input parameters. In addition, our results
predict a wide range of the orbital period (from less than one day to more
than several hundred days) for eBMSPs, which can be tested by future
observations.

###### keywords:

stars: neutron – (stars): binaries: general – (stars): pulsars: individual:
PSR J0955$-$6150

##  1 Introduction

In the standard formation theory, binary millisecond pulsars (BMSPs) originate
from low-mass X-ray binaries (LMXBs) which have undergone long time ($\sim
10^{8}-10^{10}$ yr) mass transfer (Alpar et al., 1982; Radhakrishnan &
Srinivasan, 1982; Bhattacharya & van den Heuvel, 1991). In the LMXB phase, the
neutron star (NS)’s spin is accelerated to ms periods and the binary orbit is
circularized by strong tidal interactions. As a consequence, the resultant
BMSPs are expected to reside in circular orbits (Phinney, 1992). In fact, the
majority of BMSPs with Helium white dwarf (He WD) companions in the Galactic
field have very small eccentricities in the range of $10^{-7}-10^{-3}$,
consistent with the theoretical predictions. However, recent observations
reveal the existence of a population of eccentric BMSPs (eBMSPs), which
challenges the standard formation theory. Champion et al. (2008) firstly
reported the discovery of a BMSP J1903+0327 with a high eccentricity ($e\sim
0.44$) and a G-type main-sequence companion star. The current consensus is
that this binary originated from the disruption of a triple system (Freire et
al., 2011; Portegies Zwart et al., 2011; Pijloo et al., 2012). After that,
several eBMSPs with He WD companions in the Galactic field were discovered
(see table 1 of Serylak et al., 2022, for a list of eBMSPs). The properties of
most eBMSPs can be briefly summarized as follows: (1) the eccentricities
$e\sim 0.0274-0.14$; (2) similar orbital periods $\sim 20-30$ days and WD
companion masses $\sim 0.25-0.3\,M_{\rm\odot}$; (3) fully recycled pulsars
with short spin periods $\sim 2-11$ ms; and (4) a wide range of the neutron
star (NS) masses $M_{\rm NS}\sim 1.3-1.8\,M_{\rm\odot}$.

There are five models proposed to explain (at least some of) eBMSPs in the
literature, and they can be classified into three categories, i.e., the phase
transition model (Freire & Tauris, 2014; Jiang et al., 2015), the
thermonuclear runaway burning model (Antoniadis, 2014; Han & Li, 2021), and
the resonant convection model (Ginzburg & Chiang, 2022). In the phase
transition model, the orbital eccentricity is caused by sudden mass loss due
to the phase transition of the compact object after the LMXB phase - either
accretion-induced collapse of a super-Chandrasekhar mass WD into an NS (Freire
& Tauris, 2014), or conversion of an NS into a strange quark star (Jiang et
al., 2015). In the thermonuclear runaway burning model, the eccentricity is
induced by mass loss from the proto-WDs (Antoniadis, 2014; Han & Li, 2021).
Because of unstable CNO shell burning on the surface of proto-WDs (after the
LMXB phase and prior to the WD cooling stage), part of the ejected mass may
form a circumbinary disk around the binary which dynamically interacts with
the binary and pumps up an eccentricity (Antoniadis, 2014). Han & Li (2021)
instead proposed that mass is asymmetrically lost from the proto-WDs, which
imparts a kick on the binary and makes the binary orbit to be eccentric. In
the resonant convection model, Ginzburg & Chiang (2022) argued that a coherent
resonance between the orbital period and the convective eddies in the red
giant progenitor of the WD may drive the anomalous eccentricity in eBMSPs.

Currently, all the models are subject to various uncertainties and are not
fully consistent with observations. For example, the phase transition models
predict rather specific pulsar masses, which are disfavored by the wide range
of the pulsar masses in eBMSPs. The thermonuclear runaway burning models
struggle to explain why the orbital periods of eBMSPs are concentrated in a
narrow range $\sim 20-30$ days. For the resonant convection model, it is
difficult to account for why some BMSPs with orbital periods similar to that
of eBMSPs have not undergone resonant interactions and remained in circular
orbits. Obviously, more observational clues are required to justify or put
stringent constraints on the theories.

Recent observations of PSR J0955$-$6150 carried out by Serylak et al. (2022)
showed that the binary orbit is eccentric, and more interestingly, there is
misalignment between the spin angular momentum (SAM) of the pulsar and the
orbital angular momentum (OAM). The tilt angle ($\mu$) between them is
$>4.8^{\degree}$ at 99% CI. This is the first time that a tilt angle in BMSPs
is reported. Detection of an SAM-OAM tilt angle is not surprising in other
types of double degenerate systems, such as double NSs (e.g., Weisberg et al.,
1989; Kramer, 1998), double black holes (BHs), and BH$-$NS systems (The LIGO
Scientific Collaboration et al., 2021; Abbott et al., 2021). Misalignment in
these systems is generally thought to be produced by the supernova events
during the formation of the second compact object or the dynamical process in
dense environment (e.g. Kalogera, 2000; Farr et al., 2011; Rodriguez et al.,
2016; Tauris et al., 2017; Liu & Lai, 2017; Gerosa et al., 2018). In
principle, it is possible that an SAM-OAM misalignment exists in BMSPs in
globular cluster via dynamical process. However, such misalignment is beyond
the expectation of the standard formation theory for BMSPs in the Galactic
field, because the old NSs should have experienced long duration of accretion
(Alpar et al., 1982; Bhattacharya & van den Heuvel, 1991). Though the natal
kick during the NS formation may initially excite an eccentricity and a tilt
angle, the binary orbit must have been circularized by tidal interaction and
the SAM and OAM must have been aligned by the accretion torques (Hills, 1983;
Phinney, 1992; Guillemot & Tauris, 2014). This means that the tilt angle in
PSR J0955$-$6150 was very likely produced after the mass transfer had ceased.
Unfortunately, most of the proposed formation channels of eBMSPs do not
forecast such a large tilt angle (see Serylak et al., 2022, for detailed
discussion). In addition, PSR J0955$-$6150 has a He WD with mass
$0.254(2)\,M_{\rm\odot}$, which is smaller than the theoretical prediction
according to the relation between the WD mass and orbital the period (the
$M_{\rm WD}-P_{\rm orb}$ relation, Tauris & Savonije, 1999). This is also
difficult to understand in the standard theory.

Interestingly, Serylak et al. (2022) noticed that an SAM-OAM tilt angle may be
qualitatively produced by the thermonuclear rocket model proposed by Han & Li
(2021). This model assumes that the thermonuclear flash(es) on the surface of
a proto-WD during its cooling sequence lead to asymmetrical mass ejection. As
a consequence, the proto-WD suffers mild kicks, causing the host orbit to
become eccentric and the SAM to be misaligned with the OAM. Serylak et al.
(2022) performed toy simulations of the thermonuclear rocket effect in single
flash event with the kick velocity ($v_{\rm k}$) of 8 km s-1, but failed to
obtain a large tilt angle and a long orbital period to match the observation
of PSR J0955$-$6150.

However, numerical simulations (e.g. Althaus et al., 2001; Istrate et al.,
2016) demonstrated that there are likely multiple flashes on the surface of
the proto-WDs depending on their masses. After the Roche-lobe overflow (RLO)
process, a thick hydrogen envelope (with mass $\sim 10^{-2}M_{\rm\odot}$) is
left on the outer layer of the proto-WD (Sarna et al., 2000; Althaus et al.,
2001; Benvenuto & De Vito, 2004; Istrate et al., 2016). Subsequently, the
proto-WD goes through a phase of contraction with almost constant luminosity
and increasing effective temperature. In this phase, the luminosity is mainly
supplied by the CNO burning in the shell. When the proto-WD starts cooling,
the CNO burning is not sustained because the temperature in the burning shell
is too low, and gravitational energy released due to contraction becomes the
dominant energy source. As the temperature in the hydrogen shell rises to a
critical value, the pp-chain burning is ignited, and then a flash starts with
increasing luminosity. After the flash, the proto-WD contracts with an almost
constant luminosity, and evolves following the above cycle. When the
temperature of the hydrogen shell reaches the critical value again, a second
flash occurs, and so on (See Sarna et al., 2000, for a detailed description of
the process of shell flashes). In the simulation of Istrate et al. (2016), the
maximum number of flashes can reach 26 if the effect of element diffusion
(ED)111Element diffusion is the physical mechanism for mixing of chemical
elements due to pressure, temperature and composition gradients (Thoul et al.,
1994). is considered. Thus, it is crucial to test the thermonuclear rocket
model by considering the influence of multiple flashes (or multiple kick
events).

The rest of this paper is organized as follows. In Sect. 2, we introduce the
assumptions and the simulation methods used in this work. We present our
calculated results in Sect. 3, and summarize in Sect. 4.

##  2 Assumptions and Methods

We consider binaries consisting of an NS and a proto-He WD, which is just
detached from its RL and enters the cooling sequence. The initial value of $e$
and $\mu$ are both set to be 0. The initial orbital periods ($P_{\rm orb}^{\rm
i}$) are generated from a uniform distribution. Then the WD masses are
calculated according to the $M_{\rm WD}-P_{\rm orb}$ relation with a solar
metallicity ($Z=0.02$; Tauris & Savonije, 1999). Considering the fact that the
measured masses of the NSs are in a broad range, we performed calculation with
both $M_{\rm NS}=1.4$ and $2.0\,M_{\rm\odot}$, and found that there is only
minor difference. Therefore, we only show the case of $M_{\rm
NS}=1.4\,M_{\rm\odot}$ without loss of generality.

The proto-WD is covered by a remaining hydrogen envelope. A significant part
of hydrogen is burned stably at the bottom of the envelope, while the rest is
consumed during shell flash(es) (e.g. Althaus et al., 2001; Benvenuto & De
Vito, 2004; Istrate et al., 2016). The released energy during flash(es) can
heat the outer layer material and cause the inflation of the WD’s envelope.
Han & Li (2021) argued that this energy may accelerate some outer layer
material to a velocity of $5\times 10^{3}$ km s-1, comparable to the escape
velocity of the WDs, and some of the envelope material is ejected
asymmetrically from WDs. This idea originates from Schaefer (2020) who showed
that asymmetrical mass loss during nova bursts is required to explain the
decay of the orbit after a nova burst in some cataclysmic variables.
Consequently, asymmetrical mass loss during shell flash(es) may impart a kick
to the WD and result in an eccentric orbit. By using the law of momentum
conservation, Han & Li (2021) estimated the magnitude of the kick velocity to
be around 8 km s-1, depending on the WD mass ($M_{\rm WD}$), the average
ejected mass in each flash ($\Delta M_{\rm each}$), and the velocity of the
ejecta ($V_{\rm ejct}$) (See Eq. [8]). In this work, we will examine the
influence of kicks on the orbital period, the eccentricity, and the angle
between the SAM and OAM simultaneously.

We calculate the eccentricity and the tilt angle after a flash using the
method presented in Hurley et al. (2002). Because the kick event may occur at
arbitrary position in the orbit, in each simulation we randomly choose a mean
anomaly $\mathcal{M}$ in $[0,2\pi]$ and solve the Kepler’s equation for the
eccentric anomaly $E$. The corresponding binary separation is $r=a(1-e{\rm
cos}E)$, where $a$ is the semi-major axis. It is noted that, because our
calculations involve multiple kicks and the orientation of the kick velocity
is assumed to be randomly oriented in each flash, the orientations of OAM
after several flashes are not always coplanar. Thus, the final tilt angle is
not a simple addition or subtraction of the tilt angle after each kick event.
In Appendix A.1, we present the method to calculate the tilt angle after
multiple flashes.

An important issue in the thermonuclear rocket model is the efficiency of
tidal circularization between and during flashes. The proto-WD in the inter-
flash phase is composed of a degenerated He core and a non-degenerate,
radiation-dominated hydrogen envelope with mass $\sim 10^{-2}\,M_{\rm\odot}$.
Using the MESA code (Paxton et al., 2011, 2013, 2015, 2018), we evolve a
binary with $M_{\rm NS}^{\rm i}=1.4\,M_{\rm\odot},\,M_{\rm donor}^{\rm
i}=1.4\,M_{\rm\odot}$, and $P_{\rm orb}^{\rm i}=$ 1.07 days, until the
formation of a BMSP with a 0.25 $M_{\rm\odot}$ WD companion, roughly matching
the observation of PSR J0955$-$6150\. We obtain the density and temperature
profiles from the MESA’s output files. Then we solve the equation of tidal
torque presented in Campbell (1984, Eq. 41) to evaluate the total tidal torque
$T_{\rm tot}$ and the torque exerted by the degenerate He core ($T_{\rm He,\rm
core}$). We find that $T_{\rm He,core}/T_{\rm tot}\sim 98\%$, i.e. the
envelope only contributes a tiny fraction to the total torque. Consequently,
we can safely use the tidal circularization timescale $\tau_{\rm circ}$
designed for double-degenerate binaries to estimate the efficiency of orbital
circularization in the inter-flash phase (Campbell, 1984; Hurley et al.,
2002), in which the radiative damping mechanism is adopted for the WD’s
interior. We get $\tau_{\rm circ}\approx 1.5\times 10^{26}$ yr for the above
example and $4.0\times 10^{21}$ yr for a binary consisting of a $\sim
0.20\,M_{\rm\odot}$ proto-WD. The extremely long circularization timescales
indicate that tidal circularization is negligible in the inter-flash phase,
which is consistent with the conclusion of Istrate et al. (2016).

We then consider tidal torque during the flashes. In this phase the proto-WD
can dramatically inflate and possibly fill its RL with a duration up to $\sim
10^{\rm 2}-10^{3}$ yr (e.g. Benvenuto & De Vito, 2004). The released energy
during the flashes gives rise to a pulse-driven, deep convective zone within
the hydrogen burning shell (Istrate et al., 2016). Thus, equilibrium tide is
more suitable in this situation (Yin Qin, private communication). Assume that
the WD’s spin synchronizes with the orbit, the evolution of $a$ and $e$
follows (Hut, 1981; Hurley et al., 2002),

$\displaystyle\frac{de}{dt}$ $\displaystyle=$
$\displaystyle-27\left(\frac{k}{T}\right)q\left(1+q\right)\left(\frac{R_{\rm
WD}}{a}\right)^{8}\frac{e}{\left(1-e^{2}\right)^{13/2}}$ (1)
$\displaystyle\times\left[f_{3}\left(e^{2}\right)-\frac{11}{18}\left(1-e^{2}\right)^{3/2}f_{4}\left(e^{2}\right)\right],$

and

$\displaystyle\frac{da}{dt}$ $\displaystyle=$
$\displaystyle-6\left(\frac{k}{T}\right)q\left(1+q\right)\left(\frac{R_{\rm
WD}}{a}\right)^{8}\frac{a}{\left(1-e^{2}\right)^{15/2}}$ (2)
$\displaystyle\times\left[f_{1}\left(e^{2}\right)-\left(1-e^{2}\right)^{3/2}f_{2}\left(e^{2}\right)\right],$

where $f_{n}$s are the polynomial expressions in $e^{2}$ given by Hut (1981),
$k$ the apsidal motion constant, and $T$ the timescale on which significant
changes in the orbit take place. The mass ratio $q=M_{\rm WD}/M_{\rm NS}$, and
the combined parameter of ($k/T$) can be expressed as (Rasio et al., 1996)

$\left(\frac{k}{T}\right)=\frac{2}{21}\frac{f_{\rm conv}}{{\tau}_{\rm
conv}}\frac{M_{\rm Env}}{M}\,{\rm yr}^{-1}.$ (3)

Here the convection turnover timescale

${\tau}_{\rm conv}=0.4311\left[\frac{M_{\rm Env}R_{\rm Env}\left(R_{\rm
WD}-\frac{1}{2}R_{\rm Env}\right)}{3L}\right]^{1/3}\,{\rm yr},$ (4)

where $M_{\rm Env}$ ($\sim 10^{-2}\,M_{\rm\odot}$) is the mass of the
hydrogen-rich envelope, and $R_{\rm Env}$ the depth of the convective
envelope. The luminosity $L$ of the proto-WD during the inflation phase
evolves with time, and we take $L\sim 10^{\rm 2}\,L_{\rm\odot}$ according to
the numerical calculations of Benvenuto & De Vito (2004). The numerical factor
$f_{\rm conv}\leq 1$ and is simply set to be 1 here. Note that smaller $f_{\rm
conv}$ implies longer tidal circularization time. In addition, we assume that
the proto-WD fills its RL (with radius $R_{\rm L,WD}$) during the whole flash,
and calculate the WD’s radius with Eggleton (1983)’s equation

$\frac{R_{\rm WD}}{a}=\frac{R_{\rm
L,WD}}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+{\rm ln}\left(1+q^{1/3}\right)}.$ (5)

Because the radius of the degenerate He-core is much smaller and the inflated
envelope dominates the size of the WD, we let $R_{\rm Env}=R_{\rm WD}=R_{\rm
L,WD}$. We then incorporate the above equations into our code and let the
tidal torque act on the binary for a time randomly between $10^{2}$ and
$10^{3}$ yr. Note that our treatment of $f_{\rm conv}$ and $R_{\rm WD}$ might
cause over-estimate of the tidal efficiency.

##  3 results

###  3.1 Single kick model for eBMSPs

We firstly present the predicted orbital properties of post-kick systems by
considering a single kick event, which was also investigated by Serylak et al.
(2022). Their results with a constant $v_{\rm k}=8$ km s-1 produced too small
tilt angles (always $<4.7^{\degree}$), and could not reproduce the observed
orbital period of PSR J0955$-$6150\. Here we adopt a slightly more
sophisticated model assuming that the orientation of $v_{\rm k}$ is randomly
distributed and the magnitude of $v_{\rm k}$ follows a Maxwellian distribution
with standard deviation $\sigma_{\rm k}$. We consider three cases with
$\sigma_{\rm k}=1$, 2, and 3 km s-1, but with the same ejecta mass $\Delta
M_{\rm ejct,tot}\,=10^{-3}M_{\rm\odot}$ (Han & Li, 2021), and simulate
$5\times 10^{6}$ systems with the initial WD mass in the range of
$[0.25,\,0.30]\,M_{\rm\odot}$ in each case. Fig. 1 illustrates the simulated
distributions of $P_{\rm orb}$, $e$, and $\mu$. On the top of the $P_{\rm
orb}-e$ and $P_{\rm orb}-\mu$ planes, we also show the probability
distribution functions (PDFs) of $e$ and $\mu$, respectively. The observed
BMSPs are plotted with the green stars in the top panels, including both
eBMSPs and normal BMSPs with orbital periods within $20-30$ days [that is,
PSRs J0732+2314 (Martinez et al., 2019), J1421$-$4409 (Spiewak et al., 2020),
J1709+2313 (Lewandowski et al., 2004), J1912$-$0952 (Miao et al., 2022), and
J2010+3051 (Parent et al., 2019)]. The black star represents PSR
J0955$-$6150\. The green lines in the top and bottom panels represent the
theoretical $e-P_{\rm orb}$ (Phinney, 1992) and $P_{\rm orb}-M_{\rm WD}$
(Tauris & Savonije, 1999) relations, respectively.

Figure 1: The parameter distribution of the simulated NS+He WD systems in the
single kick model. The upper, middle, and bottom rows show the distribution of
the orbital period as a function of the eccentricity, tilt angle and WD mass,
respectively. In the upper and middle panels, we also show the probability
distribution functions (PDFs) of $e$ and $\mu$, respectively. The left,
middle, and right columns represent the cases of $\sigma_{\rm k}=1$, 2 and 3
km s-1, respectively. The green stars in the upper row mark the observed
normal and eccentric BMSPs with orbital periods in the range of $\sim 20-65$
days. The black star marks PSR J0955$-$6150, and in the bottom row only eBMSPs
are plotted. The green lines in the top and bottom rows are the theoretical
$P_{\rm orb}-e$ and $P_{\rm orb}-M_{\rm WD}$ relations (Phinney, 1992; Tauris
& Savonije, 1999), respectively. The color bars on the right of the panels
represent the magnitude of $v_{\rm k}$.

The top panels in Fig. 1 demonstrate that the distribution of the eccentricity
shifts toward larger value with increasing $\sigma_{\rm k}$: the peak values
are 0.016, 0.032, and 0.056 in the cases of $\sigma_{\rm k}=$1, 2, and 3 $\rm
km\,s^{-1}$, respectively. In all the three cases the simulated results can
only partially cover the eBMSPs in the $P_{\rm orb}-e$ plane. In the middle
panels, the tilt angle also exhibits the similar tendency. Even when
$\sigma_{\rm k}=3\,\rm km\,s^{-1}$, about $3\%$ systems can reach a tilt angle
$\mu>4.8^{\rm\degree}$. Moreover, in the bottom panels, none of the simulated
systems can match the orbital period of PSR J0955$-$6150\. Therefore, we
conclude that the single kick scenario seems not to be adequate for the
formation of eBMSPs.

Figure 2: The parameter distribution of the simulated NS+He WD systems with
multiple flashes. The magnitude and orientation of the kick velocity are both
generated independently in each kick event. The kick velocity is assumed to
follow Maxwellian distribution in the upper two rows and flat distribution in
the lower two rows. For each distribution of $v_{\rm k}$, the upper and lower
panels correspond to $N=3$ and 5, respectively. The values of $v_{\rm k}$ are
limited to 5 km s-1 (colorful dots) or 10 km s-1 (grey dots). The PDFs are
shown on the top of each panel.

###  3.2 Multiple kick model for eBMSPs

As pointed out before, a proto-WD likely undergoes multiple flashes depending
on its mass. Previous numerical simulations show that the maximum number of
flashes can reach 26 if the effect of ED is considered, although in most
proto-WDs the number ($N$) of flashes is no more than $5-6$ (Istrate et al.,
2016). Thus, in this subsection we investigate the possible influence of
multiple kicks on the orbital evolution. In the following calculations, we set
$N=3$ and 5, and the total ejected mass $\Delta M_{\rm
ejct,tot}=10^{-3}\,M_{\rm\odot}$. We assume that in each flash, the ejected
mass $\Delta M_{\rm each}$ = $\Delta M_{\rm ejct,tot}/N$ and the kick velocity
is randomly orientated.

An example of our simulated results is shown in Fig. 2. Here the magnitude of
the kick velocity is assumed to follow a Maxwellian distribution (with
$\sigma_{\rm k}=1$ and 2 km s-1) or flat distribution (with $v_{\rm k}$ lying
between 0 and 10 km s-1). The upper two panels correspond to the Maxwellian
distribution of $v_{\rm k}$ with $N=3$ and 5, and the lower two panels the
flat distribution of $v_{\rm k}$ with $N=3$ and 5, respectively. The PDFs are
also demonstrated on the top of each panel. The colorful and grey dots (and
the blue and grey curves) stand for the post-kick systems with the last kick
velocity $\leq 5$ and 10 km s-1, respectively.

A comparison of Figs. 1 and 2 shows that multiple kicks indeed produce larger
eccentricities and tilt angles. For example, for Maxwellian distribution of
$v_{\rm k}$ with $\sigma_{\rm k}=2\,\rm km\,s^{-1}$ and $N=5$, 69.7% of the
systems have an eccentricity within the observed interval ($0.0274-0.1345$) of
eBMSPs, and 11.85% of the systems have a tilt angle $\mu>4.8^{\degree}$. We
also see that flat $v_{\rm k}$ distribution can produce slightly more systems
with large eccentricities and tilt angle than Maxwellian distribution, because
flat distribution tends to generate more systems with relatively high $v_{\rm
k}$. In addition, the increase in $N$ leads to larger change in the orbital
period. By comparing with the observations, we conclude that multiple flashes
with both Maxwellian and flat distributions of the kick velocity seem possible
to reproduce the properties of eBMSPs.

We then construct a number of models by changing the key parameters that
describe the kick properties, including the number of flashes ($N$) and the
kick velocity distribution $\sigma(k)$ in various flashes. The meanings of the
model names are summarized in Table 1. For example, model
$N5\\_\sigma(-2)\\_A$ means that there are 5 flashes, and in each flash, the
kick velocity independently follows flat distribution between 0 and 10 km s-1,
and model $N3\\_\sigma(+1)\\_B$ means that there are 3 flashes, and in all
flashes the kick velocity follows the same Maxwellian distribution with the
standard deviation $\sigma_{\rm k}=1$ km s-1.

In Fig. 3, we compare the simulated PDFs of the eccentricity (top panel) and
the tilt angle (bottom panel) for 19 representative models. We also list the
percentage ($f$) of systems with $e$ or $\mu$ in specific intervals confined
by the dashed lines. The green dots in the top panel represent eBMSPs and
normal BMSPs, and the red line in the bottom panel denotes the lower limit of
the tilt angle of PSR J0955$-$6150\. More detailed information about the
models and the corresponding results are presented in Table 2.

From Fig. 3 and Table 2 we see that most of the multiple kick models can more
or less cover the observed $e$ distribution of eBMSPs. However, to reach a
tilt angle $\mu>4.8\degr$ with a sufficiently high probability requires
relatively large $v_{\rm k}$ (with $\sigma_{\rm k}\geq 2$ kms-1 or $v_{\rm
k}\in[0,10]$ kms-1) and/or $N$ ($\geq 3$).

Table 1: Description of the model parameters N$i$ |  $i$ denotes the number of flashes.
---|---
$\sigma(k)$ |  if $k>0$, $v_{\rm k}$ follows Maxwellian distribution with the standard deviation $\sigma_{\rm k}=k$ $\rm km\,s^{-1}$; if $k<0$, $v_{\rm k}$ follows flat distribution, with $k=-1$ and $-2$ corresponding to $v_{\rm k}\in(0,5)$ and $(0,10)$ $\rm km\,s^{-1}$, repsectively.
Case |  A: the $v_{\rm k}$ distribution is independently set in each flash. B: the $v_{\rm k}$ distributions in all flashes take the same form.

Figure 3: The PDFs of $e$ (upper panel) and $\mu$ (lower panel) for 19
representative models. In the upper panel, the green stars represent the
observed eccentricities of BMSPs. In the lower panel, the red line represents
the lower limit of the tilt angle of PSR J0955$-$6150\. The values of $f$
denote the percentages of the systems located in specific intervals.

Figure 4: Similar to Fig. 3, but with the asymmetry parameter $f_{\rm a}$
considered in the left 10 columns, and with fitted relations considered in the
right 9 columns. The inset in the lower panel shows the PDF in model
$N5\\_\sigma(-2)\\_B+fa\\_U$.

###  3.3 Multiple kick model for both BMSPs and eBMSPs

A successful kick model should be able to simultaneously account for both
eBMSPs and normal BMSPs. While the above models seem to be able to explain the
properties of eBMSPs, none of them can cover the eccentricity distribution of
normal BMSPs. This means that the realistic situations of shell flashes are
more complicated than our assumptions about the kick velocity, that is, the
kick velocity may depend on not only the ejecta mass but also the asymmetry of
mass ejection. So we introduce an asymmetry parameter $f_{\rm a}$ which
quantifies the extent of asymmetry in the mass loss processes during shell
flashes. The case of $f_{\rm a}=1$ corresponds to fully asymmetric ejection,
while $f_{\rm a}=0$ represents fully symmetric mass ejection without any kick.
Then the modified kick velocity is $f_{\rm a}v_{\rm k}$. Since the form of
$f_{\rm a}$ is unclear, we tentatively consider two possibilities: (1)
$\log\,f_{\rm a}$ follows a flat distribution in the range of $[-5,0]$, and
(2) $\log f_{\rm a}=C(1-x)^{3}$, where $C=-5$ and $x$ is a random number in
$[0,1]$. While in the former case the $f_{\rm a}$ distribution is dominated by
relatively small values, in the latter case it is clustered around both small
and large values. We add “fa” and “fa_U” to the original model name to
identify them, respectively.

We show the calculated results for 10 representative models with asymmetrical
mass ejection considered in the left 10 columns of Fig. 4 and in Table B1.
Comparing with Fig. 3 we find that the distributions of the eccentricity cover
a wider range than predicted by models without $f_{\rm a}$ considered, but the
tilt angle is clustered smaller value. Among all the models, model
$N5\\_\sigma(-2)\\_B+fa\\_U$, seems to simultaneously match the eccentricities
of eBMSPs and normal BMSPs and the tilt angle of PSR J0955$-$6150\. The inset
in the lower panel shows the PDF of $\mu$ in detail.

###  3.4 The effect of flash number and ejected mass

We finally discuss the possible dependence of the flash number and ejected
mass on the WD mass, rather taking them as fixed values. Depending on
metallicity and ED, different WDs may experience different times of flashes
(Althaus et al., 2001; Istrate et al., 2016). Generally $N$ is smaller for
more massive WDs with the same metallicity (Istrate et al., 2016). In
addition, the ejected mass during flashes may be different for different WDs.
For example, Benvenuto & De Vito (2004) found that only $\sim$10% of the
hydrogen envelope mass ($M_{\rm Env,det}$) is lost for a $0.21\,M_{\rm\odot}$
WD during the flashes, while $\sim 43\%$ of $M_{\rm Env,det}$ is lost for a
$0.3\,M_{\rm\odot}$ WD in Iben & Tutukov (1986).

Basing on the calculated data provided in Table A1 of Istrate et al. (2016),
we derive an empirical $N-M_{\rm WD}$ relation by utilizing polynomial fitting
(See Appendix A.2 for details). Then the average lost mass $\Delta M_{\rm
each}$ during each flash is taken to be

$\Delta M_{\rm each}=\frac{\Delta M_{\rm ejct,tot}}{N(M_{\rm WD})},$ (6)

where $\Delta M_{\rm ejct,tot}$ is also fitted as a function of $M_{\rm WD}$
by using the data of $M_{\rm Env,det}$ in Istrate et al. (2016), i.e.,

$\Delta M_{\rm ejct,tot}(M_{\rm WD})=f_{\rm ejct}M_{\rm Env,det}(M_{\rm WD}).$
(7)

Here $f_{\rm ejct}$ denotes the percentage of mass loss to $M_{\rm Env,det}$
during flashes, and is set to be randomly in the range of $[10\%,50\%]$
according to Iben & Tutukov (1986) and Benvenuto & De Vito (2004). Then the
kick velocity is calculated with (see also Han & Li, 2021)

$v_{\rm k}=8\,{\rm km\,s^{-1}}f_{\rm a}\left(\frac{\Delta M_{\rm
each}}{10^{-3}M_{\sun}}\right)\left(\frac{V_{\rm ejct}}{\rm 2.5\times
10^{3}\,km\,s^{-1}}\right)\left(\frac{M_{\rm WD}}{0.3M_{\sun}}\right)^{-1},$
(8)

where the ejecta velocity $V_{\rm ejct}$ is set to be $2.5\times 10^{3}$ km
s-1. We assume that the magnitude of $v_{\rm k}$ is the same for all flashes
but its orientation is randomly set during each flash. We consider three
cases: (1) the metallicity $Z=0.02$ without ED (named as Z02_ED0), (2)
$Z=0.02$ with ED (Z02_ED1), and (3) $Z=0.001$ without ED (Z001_ED0). A log-
normal distribution is adopted for $f_{\rm a}$. The simulated results are
presented in the right 9 columns of Fig. 4 and Table B2. In the model names,
$f10$, $f30$, and $f50$ correspond to $f_{\rm ejct}=10\%$, $30\%$, and $50\%$,
respectively. Note that here we include all systems with the occurrence of
shell flash(es). Therefore, the range of $M_{\rm WD}$ is significantly wider
than the adopted range of $[0.25,0.30]\,M_{\sun}$ in previous discussion.
Specifically, $M_{\rm WD}$ is in $[0.213,0.305]\,M_{\sun}$,
$[0.167,0.393]\,M_{\sun}$ and $[0.252,0.350]\,M_{\sun}$ for models of Z02_ED0,
Z02_ED1, Z001_ED0, respectively. Thus, we plot all BMSPs with known
eccentricities for comparison.

The right 9 columns in Fig. 4 show the following features.

(1) When $f_{\rm ejct}=10\%$, the kick velocities can barely explain the
observed eccentricities of eBMSPs. The predicted tilt angles ($<3^{\degree}$)
are also smaller than the lower limit of $\mu$ for PSR J0955$-$6150\.
Increasing $f_{\rm ejct}$ can produce more systems with eccentric orbits and
larger tilt angles. We find that models Z02_ED0_f50_B+fa, Z001_ED0_f30_B+fa,
and Z001_ED0_f50_B+fa seem to better match the observational properties of
both eBMSPs and normal BMSPs.

(2) When ED is considered, both the eccentricities and tilt angles shift to
smaller values compared with those without ED, due to relatively smaller
ejected mass and kick velocity.

(3) Models with lower metallicity generally produce larger eccentricities and
tilt angles under the same $f_{\rm ejct}$. This difference originates from the
more massive envelope (i.e., larger $M_{\rm Env,det}$) for lower metallicity
stars (Istrate et al., 2016), which in turn results in larger kick velocity
(See Fig. 7).

##  4 Summary

In this work, we have investigated the feasibility of the thermonuclear rocket
model for the formation of eBMSPs. Here, we focus on how to explain the tilt
angel between the SAM and the OAM and the orbital period (or the WD mass)
deviation from the theoretical prediction for PSR J0955$-$6150\. We also
explore the possibility to simultaneously explain the overall eccentricity
distribution of both eBMSPs and normal BMSPs under the thermonuclear rocket
model.

We first let all WDs have the same number $N$ of flashes and the same ejecta
mass $\Delta M_{\rm ejec,tot}$ during flashes, and artificially limit the mass
range of WDs within that of observed eBMSPs. Under these simplified
assumptions, we show that the thermonuclear rocket model with multiple kicks
could better reproduce the properties of both eBMSPs than with single kick.
However, to account for the normal BMSPs simultaneously, one needs to
carefully take into consideration the effect of asymmetry of mass ejection.

We then discuss the more realistic case based on the numerical data of shell
flashes on proto-WDs in previous studies. In this case, both $N$ and $\Delta
M_{\rm ejec,tot}$ depend on $M_{\rm WD}$, metallicity and ED. We fit the
relations among these factors and incorporate them into our simulations. We
find that under such situation, a few models with relatively large (30% - 50%)
mass ejection fraction can potentially match the observations. Moreover, our
results predict a wide range of the orbital period (from less than one day to
more than several hundred days) for eBMSPs, and some of them may have a tilt
angle up to $\sim 10\degr$. We expect that these predictions to be tested by
future observations.

## Acknowledgments

We thank the referee, Pablo Marchant, for constructive comments and
suggestions that have helped improve the manuscript. This work was supported
by the National Key Research and Development Program of China
(2021YFA0718500), the Natural Science Foundation of China under grant No.
12041301, 12121003, and Project U1838201 supported by NSFC and CAS.

## Data Availability

All data underlying this article will be shared on reasonable request to the
corresponding authors.

## References

  * Alpar et al. (1982) Alpar, M. A., Cheng, A. F., Ruderman, M. A., et al. 1982, Nature, 300, 728. doi:10.1038/300728a0
  * Radhakrishnan & Srinivasan (1982) Radhakrishnan, V. & Srinivasan, G. 1982, Current Science, 51, 1096
  * Bhattacharya & van den Heuvel (1991) Bhattacharya, D. & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1. doi:10.1016/0370-1573(91)90064-S
  * Phinney (1992) Phinney, E. S. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 39. doi:10.1098/rsta.1992.0084
  * Champion et al. (2008) Champion, D. J., Ransom, S. M., Lazarus, P., et al. 2008, Science, 320, 1309. doi:10.1126/science.1157580
  * Freire et al. (2011) Freire, P. C. C., Bassa, C. G., Wex, N., et al. 2011, MNRAS, 412, 2763. doi:10.1111/j.1365-2966.2010.18109.x
  * Portegies Zwart et al. (2011) Portegies Zwart, S., van den Heuvel, E. P. J., van Leeuwen, J., et al. 2011, ApJ, 734, 55. doi:10.1088/0004-637X/734/1/55
  * Pijloo et al. (2012) Pijloo, J. T., Caputo, D. P., & Portegies Zwart, S. F. 2012, MNRAS, 424, 2914. doi:10.1111/j.1365-2966.2012.21431.x
  * Deneva et al. (2013) Deneva, J. S., Stovall, K., McLaughlin, M. A., et al. 2013, ApJ, 775, 51. doi:10.1088/0004-637X/775/1/51
  * Serylak et al. (2022) Serylak, M., Venkatraman Krishnan, V., Freire, P. C. C., et al. 2022, arXiv:2203.00607
  * Serylak et al. (2022) Serylak, M., Venkatraman Krishnan, V., Freire, P. C. C., et al. 2022, A&A, 665, A53. doi:10.1051/0004-6361/202142670
  * Freire & Tauris (2014) Freire, P. C. C. & Tauris, T. M. 2014, MNRAS, 438, L86. doi:10.1093/mnrasl/slt164
  * Jiang et al. (2015) Jiang, L., Li, X.-D., Dey, J., et al. 2015, ApJ, 807, 41. doi:10.1088/0004-637X/807/1/41
  * Antoniadis (2014) Antoniadis, J. 2014, ApJ, 797, L24. doi:10.1088/2041-8205/797/2/L24
  * Han & Li (2021) Han, Q. & Li, X.-D. 2021, ApJ, 909, 161. doi:10.3847/1538-4357/abdd21
  * Ginzburg & Chiang (2022) Ginzburg, S. & Chiang, E. 2022, MNRAS, 509, L1. doi:10.1093/mnrasl/slab110
  * Kalogera (2000) Kalogera, V. 2000, ApJ, 541, 319. doi:10.1086/309400
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170. doi:10.3847/1538-4357/aa7e89
  * Althaus et al. (2001) Althaus, L. G., Serenelli, A. M., & Benvenuto, O. G. 2001, MNRAS, 323, 471. doi:10.1046/j.1365-8711.2001.04227.x
  * Istrate et al. (2016) Istrate, A. G., Marchant, P., Tauris, T. M., et al. 2016, A&A, 595, A35. doi:10.1051/0004-6361/201628874
  * Benvenuto & De Vito (2004) Benvenuto, O. G. & De Vito, M. A. 2004, MNRAS, 352, 249. doi:10.1111/j.1365-2966.2004.07918.x
  * Sarna et al. (2000) Sarna, M. J., Ergma, E., & Gerškevitš-Antipova, J. 2000, MNRAS, 316, 84. doi:10.1046/j.1365-8711.2000.03503.x
  * Schaefer (2020) Schaefer, B. E. 2020, MNRAS, 492, 3343. doi:10.1093/mnras/stz3424
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897. doi:10.1046/j.1365-8711.2002.05038.x
  * Hills (1983) Hills, J. G. 1983, ApJ, 267, 322. doi:10.1086/160871
  * Tauris & Savonije (1999) Tauris, T. M. & Savonije, G. J. 1999, A&A, 350, 928
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3. doi:10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4. doi:10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15. doi:10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34. doi:10.3847/1538-4365/aaa5a8
  * Zahn (1975) Zahn, J.-P. 1975, A&A, 41, 329
  * Lorimer et al. (2021) Lorimer, D. R., Kawash, A. M., Freire, P. C. C., et al. 2021, MNRAS, 507, 5303. doi:10.1093/mnras/stab2474
  * Guillemot & Tauris (2014) Guillemot, L. & Tauris, T. M. 2014, MNRAS, 439, 2033. doi:10.1093/mnras/stu082
  * Farr et al. (2011) Farr, W. M., Kremer, K., Lyutikov, M., et al. 2011, ApJ, 742, 81. doi:10.1088/0004-637X/742/2/81
  * Rodriguez et al. (2016) Rodriguez, C. L., Zevin, M., Pankow, C., et al. 2016, ApJ, 832, L2. doi:10.3847/2041-8205/832/1/L2
  * Liu & Lai (2017) Liu, B. & Lai, D. 2017, ApJ, 846, L11. doi:10.3847/2041-8213/aa8727
  * Camilo et al. (2015) Camilo, F., Kerr, M., Ray, P. S., et al. 2015, ApJ, 810, 85. doi:10.1088/0004-637X/810/2/85
  * Octau et al. (2018) Octau, F., Cognard, I., Guillemot, L., et al. 2018, A&A, 612, A78. doi:10.1051/0004-6361/201731928
  * Edwards & Bailes (2001) Edwards, R. T. & Bailes, M. 2001, ApJ, 553, 801. doi:10.1086/320986
  * Barr et al. (2013) Barr, E. D., Champion, D. J., Kramer, M., et al. 2013, MNRAS, 435, 2234. doi:10.1093/mnras/stt1440
  * Barr et al. (2017) Barr, E. D., Freire, P. C. C., Kramer, M., et al. 2017, MNRAS, 465, 1711. doi:10.1093/mnras/stw2947
  * Knispel et al. (2015) Knispel, B., Lyne, A. G., Stappers, B. W., et al. 2015, ApJ, 806, 140. doi:10.1088/0004-637X/806/1/140
  * Zhu et al. (2019) Zhu, W. W., Freire, P. C. C., Knispel, B., et al. 2019, ApJ, 881, 165. doi:10.3847/1538-4357/ab2bef
  * Stovall et al. (2019) Stovall, K., Freire, P. C. C., Antoniadis, J., et al. 2019, ApJ, 870, 74. doi:10.3847/1538-4357/aaf37d
  * Martinez et al. (2019) Martinez, J. G., Gentile, P., Freire, P. C. C., et al. 2019, ApJ, 881, 166. doi:10.3847/1538-4357/ab2877
  * Spiewak et al. (2020) Spiewak, R., Flynn, C., Johnston, S., et al. 2020, MNRAS, 496, 4836. doi:10.1093/mnras/staa1869
  * Lewandowski et al. (2004) Lewandowski, W., Wolszczan, A., Feiler, G., et al. 2004, ApJ, 600, 905. doi:10.1086/379923
  * Parent et al. (2019) Parent, E., Kaspi, V. M., Ransom, S. M., et al. 2019, ApJ, 886, 148. doi:10.3847/1538-4357/ab4f85
  * Miao et al. (2022) Miao, C. C., Zhu, W. W., Li, D., et al. 2022, arXiv:2205.04232
  * Antoniadis et al. (2016) Antoniadis, J., Kaplan, D. L., Stovall, K., et al. 2016, ApJ, 830, 36. doi:10.3847/0004-637X/830/1/36
  * Weisberg et al. (1989) Weisberg, J. M., Romani, R. W., & Taylor, J. H. 1989, ApJ, 347, 1030. doi:10.1086/168193
  * Kramer (1998) Kramer, M. 1998, ApJ, 509, 856. doi:10.1086/306535
  * Gerosa et al. (2018) Gerosa, D., Berti, E., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 98, 084036. doi:10.1103/PhysRevD.98.084036
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv:2111.03634
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 915, L5. doi:10.3847/2041-8213/ac082e
  * Driebe et al. (1999) Driebe, T., Blöcker, T., Schönberner, D., et al. 1999, A&A, 350, 89
  * Thoul et al. (1994) Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828. doi:10.1086/173695
  * Iben & Tutukov (1986) Iben, I. & Tutukov, A. V. 1986, ApJ, 311, 742. doi:10.1086/164812
  * Campbell (1984) Campbell, C. G. 1984, MNRAS, 207, 433. doi:10.1093/mnras/207.3.433
  * Tout et al. (1997) Tout, C. A., Aarseth, S. J., Pols, O. R., et al. 1997, MNRAS, 291, 732. doi:10.1093/mnras/291.4.732
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., et al. 2005, AJ, 129, 1993. doi:10.1086/428488
  * Rasio et al. (1996) Rasio, F. A., Tout, C. A., Lubow, S. H., et al. 1996, ApJ, 470, 1187. doi:10.1086/177941
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368. doi:10.1086/160960

## Appendix A

###  A.1 The tilt angle after multiple kicks

The geometry of the binary is shown in Fig. 5. The NS is at the coordinate
origin, and the orbit is on the $xy$ plane. We assume that the unit vector
$\hat{s}$ of the spin momentum of the NS is initially aligned with the orbital
angular momentum $\hat{h}_{0}$ and the $z$-axis, i.e.,
$\hat{s}=\hat{h}_{0}=(0,0,1)$.

After the first kick, the angle $\mu$ between $\hat{s}$ and $\hat{h}_{1}$ is
$\mu_{1}=\theta_{1}=\langle\hat{s},\hat{h}_{1}\rangle$. The toroidal angle of
$\hat{h}_{1}$ is arbitrarily set, and we assume that ${\hat{h}}_{1}$ is in the
$yOz$ plane.

After the second kick, the orbital angular momentum vector becomes
$\hat{h}_{2}$. Using the method in Appendix A1 of Hurley et al. (2002), we can
obtain the angle between $\hat{h}_{1}$ and $\hat{h}_{2}$ , i.e., $\theta_{2}$.
Since the angular momentum vector $\hat{h}_{2}$ is not coplanar with
$\hat{h}_{0}$ and $\hat{h}_{1}$, the final tilt angle is not linear addition
or subtraction between $\theta_{1}$ and $\theta_{2}$. We derive $\mu$ after
the second kick in the following way. We first rotate $\hat{h}_{1}$ by the
angle $\theta_{2}$ to get $\hat{h}^{\prime}_{2}$ in the $yOz$ plane. Then, we
use the Rodrigues’ Rotation Formula

$\hat{h}_{2}=\cos\varphi\hat{h}^{\prime}_{2}+(1-\cos\varphi)(\hat{h}^{\prime}_{2}\cdot\hat{h}_{1})\hat{h}_{1}+\hat{h}_{1}\times\hat{h}^{\prime}_{2}\sin\varphi,$
(9)

to rotate $\hat{h}^{\prime}_{2}$ along $\hat{h}_{1}$ by an arbitrary angle
from a uniform distribution of $[0,2\pi]$ and finally get $\hat{h}_{2}$. So,
the angle between the current orbital angular momentum and the spin angular
momentum of the NS is $\mu=\langle\hat{s},\hat{h}_{2}\rangle$ (which is
$\leq\theta_{1}+\theta_{2}$). For multiple kicks, $\mu$ can be obtained by
repeating the above steps.

###  A.2 Fitting the $N-M_{\rm WD}$ and $\Delta M_{\rm ejct,tot}-M_{\rm WD}$
relations

Based on the data provided in Table A1 of Istrate et al. (2016), we utilize
the polynomial function

$F(x)=\sum_{\rm j=0}^{\rm n}w_{\rm j}x^{\rm j}$ (10)

to fit the $N-M_{\rm WD}$ and $\Delta M_{\rm ejct,tot}-M_{\rm WD}$ relations
with the polyfit function in the python code, where $n=5$ for $N(M_{\rm WD})$
and 8 for $M_{\rm ejct,tot}(M_{\rm WD})$. An example is demonstrated in Fig.
6.

Fig. 7 shows several examples of the $v_{\rm k}$ distribution by use of the
empirical $N-M_{\rm WD}$ and $\Delta M_{\rm ejct,tot}-M_{\rm WD}$ relations
and Eq. (8). Here we adopt a log-normal $f_{\rm a}$ distribution and $f_{\rm
ejct}=30\%$ (left) and 50% (right).

Figure 5: The geometry of the model for calculating $\mu$ after multiple
kicks. Figure 6: An example of the empirical $N-M_{\rm WD}$ (left) and $\Delta
M_{\rm ejct,tot}-M_{\rm WD}$ (right) relations for model of Z02_ED1.  Figure
7: The obtained $v_{\rm k}$ distribution by using the fitting result of
Appendix. A.2 and Eq.( 8) to calculate $v_{\rm k}$, where a log-normal $f_{\rm
a}$ distribution has been taken into account. The left and right panels show
the case with $f_{\rm ejct}=$ 30% and 50%, respectively. In each panel, the
green, black and red lines correspond to the models of Z001_ED0, Z02_ED0,
Z02_ED1, respectively.

## Appendix B

Table 2 presents a list of the model parameters and the calculated results.

Table 2: Properties and predictions of multiple kick Models.

Name Properties |  $v_{\rm k}$ distribution |  $v_{\rm k}$ limit(1) | $N$ |  $v_{\rm k}$ in flashes | $f_{\rm a}^{(2)}$ | $e\leq 10^{-3}$ | $e\geq 10^{-2}$ | $e\geq 10^{-1}$ | $\mu\geq 4.8^{\degree}$ | $\mu\geq 10^{\degree}$
---|---|---|---|---|---|---|---|---|---|---
N1_$\sigma(+1)$_A | Maxwell: $\sigma_{\rm k}=1^{(1)}$ | … | 1 |  | No | 0.2% | 84.0% |  |  |
N1_$\sigma(+2)$_A | Maxwell: $\sigma_{\rm k}$ =2 | … | 1 |  | No |  | 95.7% | 5.2% | 0.1% |
N1_$\sigma(+3)$_A | Maxwell: $\sigma_{\rm k}$ =3 | … | 1 |  | No |  | 98.0% | 21.1% | 2.9% |
N3_$\sigma(+1)$_A | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 3 | independent | No | 0.1% | 94.7% | 1.6% |  |
N3_$\sigma(+2)$_A | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 3 | independent | No |  | 98.6% | 29.0% | 4.1% |
N5_$\sigma(+1)$_A | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 5 | independent | No |  | 96.8% | 5.9% | 0.1% |
N5_$\sigma(+2)$_A | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 5 | independent | No |  | 99.2% | 44.8% | 111.9% | 0.1%
N3_$\sigma(-1)$_A | Flat | 5 | 3 | independent | No |  | 97.8% | 21.2% | 0.4% |
N3_$\sigma(-2)$_A | Flat | 10 | 3 | independent | No |  | 99.4% | 62.3% | 18.1% | 0.3%
N5_$\sigma(-1)$_A | Flat | 5 | 5 | independent | No |  | 98.6% | 34.9% | 2.4% |
N5_$\sigma(-2)$_A | Flat | 10 | 5 | independent | No |  | 99.6% | 72.3% | 34.3% | 2.3%
N3_$\sigma(+1)$_B | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 3 | unchanged | No | 0.5% | 81.7% | 2.2% |  |
N3_$\sigma(+2)$_B | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 3 | unchanged | No | 0.1% | 93.3% | 18.3% | 1.4% |
N5_$\sigma(+1)$_B | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 5 | unchanged | No | 0.4% | 85.5% | 5.5% | 0.1% |
N5_$\sigma(+2)$_B | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 5 | unchanged | No | 0.1% | 94.7% | 26.9% | 4.4% | 0.1%
N3_$\sigma(-1)$_B | Flat | 5 | 3 | unchanged | No | 1.2% | 87.5% | 20.3% | 1.4% |
N3_$\sigma(-2)$_B | Flat | 10 | 3 | unchanged | No | 0.6% | 93.7% | 47.5% | 18.8% | 1.3%
N5_$\sigma(-1)$_B | Flat | 5 | 5 | unchanged | No | 1.1% | 89.4% | 28.5% | 5.1% |
N5_$\sigma(-2)$_B | Flat | 10 | 5 | unchanged | No | 0.5% | 94.6% | 54.5% | 29.3% | 5.2%
N3_$\sigma(+1)$_B+fa | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 3 | unchanged | Yes | 65.8% | 9.8% |  |  |
N3_$\sigma(+2)$_B+fa | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 3 | unchanged | Yes | 58.2% | 16.5% | 0.7% |  |
N5_$\sigma(+1)$_B+fa | Maxwell: $\sigma_{\rm k}$ =1 | 5 | 5 | unchanged | Yes | 63.9% | 11.7% | 0.1% |  |
N5_$\sigma(+2)$_B+fa | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 5 | unchanged | Yes | 56.4% | 18.6% | 1.2% | 0.1% |
N3_$\sigma(-1)$_B+fa | Flat | 5 | 3 | unchanged | Yes | 59.6% | 15.9% | 0.8% |  |
N3_$\sigma(-2)$_B+fa | Flat | 10 | 3 | unchanged | Yes | 52.2% | 22.8% | 3.3% | 0.6% |
N5_$\sigma(-1)$_B+fa | Flat | 5 | 5 | unchanged | Yes | 57.8% | 17.8% | 1.4% | 0.1% |
N5_$\sigma(-2)$_B+fa | Flat | 10 | 5 | unchanged | Yes | 50.3% | 24.9% | 4.5% | 1.3% | 0.1%
N5_$\sigma(+2)$_B+fa_U | Maxwell: $\sigma_{\rm k}$ =2 | 10 | 5 | unchanged | “U”-like | 31.8% | 51.6 % | 19.6% | 8.4% | 1.1%
N5_$\sigma(-2)$_B+fa_U | Flat | 10 | 5 | unchanged | “U”-like | 34.3% | 46.7% | 7.8% | 1.0% |

  * 1

in the unit of $\rm km\,s^{-1}$.

  * 2

$f_{\rm a}$ follows a log-normal distribution unless otherwise specified.

Table 3: Properties and predictions of multiple kick Models based on the fitted $N-M_{\rm WD}$ and $\Delta M_{\rm ejct,tot}-M_{\rm WD}$ relations.  Name Properties |  $v_{\rm k}$ and $N$ | $Z$ |  element diffusion | $f_{\rm eject}$ |  $v_{\rm k}$ in flashes |  $f_{\rm a}$ | $e\leq 10^{-3}$ | $e\geq 10^{-2}$ | $e\geq 10^{-1}$ | $\mu\geq 4.8^{\degree}$ | $\mu\geq 10^{\degree}$
---|---|---|---|---|---|---|---|---|---|---|---
$\rm Z02\\_ED0\\_f10\\_B+fa$ | Fitting | 0.02 | No | 10% | unchanged |  Yes | 59.7% | 15.1% |  |  |
$\rm Z02\\_ED0\\_f30\\_B+fa$ | Fitting | 0.02 | No | 30% | unchanged |  Yes | 44.1% | 27.1% | 2.8% | 0.6% |
$\rm Z02\\_ED0\\_f50\\_B+fa$ | Fitting | 0.02 | No | 50% | unchanged |  Yes | 16.6% | 32.7% | 8.0% | 5.3% |
$\rm Z02\\_ED1\\_f10\\_B+fa$ | Fitting | 0.02 | Yes | 10% | unchanged |  Yes | 61.8% | 13.6% |  |  |
$\rm Z02\\_ED1\\_f30\\_B+fa$ | Fitting | 0.02 | Yes | 30% | unchanged |  Yes | 48.8% | 25.6% | 4.2% | 0.3% |
$\rm Z02\\_ED1\\_f50\\_B+fa$ | Fitting | 0.02 | Yes | 50% | unchanged |  Yes | 38.6% | 31.2% | 9.3% | 2.2% | 0.1%
$\rm Z001\\_ED0\\_f10\\_B+fa$ | Fitting | 0.001 | No | 10% | unchanged |  Yes | 54.2% | 20.4% | 0.1% |  |
$\rm Z001\\_ED0\\_f30\\_B+fa$ | Fitting | 0.001 | No | 30% | unchanged |  Yes | 24.5% | 32.4% | 8.4% | 4.5% | 0.1%
$\rm Z001\\_ED0\\_f50\\_B+fa$ | Fitting | 0.001 | No | 50% | unchanged |  Yes | 7.8% | 38.0% | 14.0% | 9.3% | 2.1%


# MedConceptsQA: Open Source Medical Concepts QA Benchmark

Ofir Ben Shoham
Department of Software and Information Systems Engineering
Ben-Gurion University of the Negev
<EMAIL_ADDRESS>
&Nadav Rappoport
Department of Software and Information Systems Engineering
Ben-Gurion University of the Negev
<EMAIL_ADDRESS>

###### Abstract

We present MedConceptsQA, a dedicated open source benchmark for medical
concepts question answering. The benchmark comprises of questions of various
medical concepts across different vocabularies: diagnoses, procedures, and
drugs. The questions are categorized into three levels of difficulty: easy,
medium, and hard. We conducted evaluations of the benchmark using various
Large Language Models. Our findings show that pre-trained clinical Large
Language Models achieved accuracy levels close to random guessing on this
benchmark, despite being pre-trained on medical data. However, GPT-4 achieves
an absolute average improvement of nearly 27%-37% (27% for zero-shot learning
and 37% for few-shot learning) when compared to clinical Large Language
Models. Our benchmark serves as a valuable resource for evaluating the
understanding and reasoning of medical concepts by Large Language Models. Our
benchmark is available at
https://huggingface.co/datasets/ofir408/MedConceptsQA

##  1 Introduction

Large Language Models (LLMs) are trained on huge datasets and comprising
billions of parameters. LLMs have demonstrated remarkable effectiveness across
a range of language-related tasks, including question-answering and text
generation [1]. Their applicability extends beyond traditional language
domains to areas such as robotics, cybersecurity, law, and also medicine [2].

Clinical LLMs (CLLMs) are LLMs that are trained on medical datasets. To name a
few: Meditron [3], BioMistral [4], BioMedGPT [5], GatorTron [6], and Meerkat
[7]. CLLMs have shown their effectiveness to be used in tasks such as clinical
text classification [8], medical chatbots, healthcare education, and clinical
text generation [9].

Clinical data is frequently represented using standardized medical codes
rather than natural language descriptions. Therefore, it’s essential that
CLLMs understand the meaning of these medical codes and their differences. One
way to test the understanding of CLLMs is to evaluate their ability to
interpret clinical codes.

Existing clinical benchmarks are available for evaluating CLLMs. For example,
BioASQ-QA [10] is a manually curated corpus of Biomedical Question Answering
(QA) for documented retrieval, text snippets extraction, and summarization.
EmrQA [11] is a corpus for QA on Electronic Medical Records based on clinical
notes. PubMedQA [12] is a dataset that contains biomedical QA based on PubMed
abstracts. In addition to QA datasets, there are clinical datasets in other
fields such as classification [13] and summarization [14, 15]. Soroush et al.
[16] used real clinical visit data to provide LLMs with a code description and
prompted them to generate a billing code. They found that all general-purpose
LLMs tested performed poorly at this task. Notably, they did not evaluate
CLLMs that were trained on medical data and they only used zero-shot. We used
a different approach for evaluation: our proposed benchmark includes a broader
set of questions and answers relating to medical concepts and also the
differences between them, spanning various levels of difficulty. Furthermore,
our evaluation encompasses not only general-purpose LLMs but also CLLMs
specialized for the medical domain. Unlike their focus on a single hospital
site, our medical coding evaluation covers a wider scope and also contains
drugs medical concepts. Additionally, our evaluation also incorporates few-
shot learning learning and not only zero-shot learning. Moreover, Soroush et
al. generated code from natural language descriptions (a task that is
currently performed by human coders), while our study evaluates the "clinical
understanding" of LLMs by investigate the inverse direction: selection the
right description given a clinical code. In contrast to their approach, which
assessed the capability of models to potentially replace human coders, our
goal is to provide a more general estimation of the clinical reasoning
abilities of LLMs.

In this study, we introduce MedConceptsQA, a benchmark for evaluating
understanding and reasoning abilities in the medical domain. MedConceptsQA
comprises over 800,000 questions and answers covering medical concepts,
including ICD10 and ICD9 diagnoses codes, ICD9-PROC and ICD10-PROC procedures
codes, and ATC drug codes. The benchmark includes questions categorized into
three difficulty levels: easy, medium, and hard.

In addition, we evaluate different LLMs on the MedConceptsQA benchmark using
zero-shot and few-shot learning. Surprisingly, we find that the current state-
of-the-art CLLMs, as well as general LLMs except two, achieve performance
levels comparable to Random guessing. Even medical fine-tuned models were not
better than Random guessing. However, GPT-4 [17] outperforms all medical LLMs,
despite its primary focus being a general LLM rather than specifically
tailored to the medical domain like the others.

Contributions: In summary, we make three main contributions:

  1. 1.

We propose a challenging open-source benchmark for evaluating LLMs on their
understanding and reasoning of medical concepts. Additionally, our evaluation
code is also open-source and available for any use.

  2. 2.

Due to our large number of examples in the benchmark (more than 800,000), our
benchmark can be utilized for training LLMs to comprehend medical concepts and
the distinctions between them, making it suitable for techniques such as
instruction-tuning [18].

  3. 3.

We demonstrate that GPT-3.5 and GPT-4 outperform the current state-of-the-art
CLLMs in understanding and reasoning about medical concepts, despite CLLMs
being specifically designed for the medical domain.

##  2 MedConceptsQA Benchmark

The MedConceptsQA benchmark consists of multiple question-answering medical
concepts datasets, categorized into three vocabularies: Diagnoses, Procedures,
and Drugs. Within the Diagnoses category, it includes questions related to
ICD9-CM and ICD10-CM medical codes. For Procedures, it covers ICD9-PROC and
ICD10-PROC codes, and ATC for Drugs. Each vocabulary is further divided into
three levels of difficulty: easy, medium, and hard. Table 1 illustrates the
distribution of questions within the benchmark.

| Vocab | Level | Questions num
---|---|---|---
| ATC | easy | 6440
medium | 6440
hard | 5938
| ICD10-CM | easy | 94580
medium | 81757
hard | 88013
| ICD10-PROC | easy | 190987
medium | 190987
hard | 88582
| ICD9-CM | easy | 17736
medium | 17736
hard | 16858
| ICD9-PROC | easy | 4670
medium | 4670
hard | 4438
Table 1: The number of generated questions across different vocabularies and
difficulty levels in the benchmark.

Each question is about a single medical code and contains four optional
answers. The options are four descriptions of medical codes, where only one is
the true description of the given medical codes, and the others are randomly
chosen according to the category of difficulty of the specific level. The
order of the options and the placement of correct answers are chosen randomly.

###  2.1 Difficulty Levels Creation

We represent the medical code vocabulary hierarchy as an undirected graph
using PyHealth [19]. Each vocabulary is divided into three difficulty levels:
easy, medium, and hard. The distinction lies in the options provided for the
correct answer. For the easy level, the options are randomly chosen from all
the medical concept codes within the vocabulary. For the medium level, we
select alternative medical codes within a distance of three, four, or five
edges from each node in the vocabulary. For the hard level, the required
distance is reduced to two edges. This means that for the hard level, the
candidates are closely related, sharing a common parent or a node and its
grand-parent or a node and its grand-child. For example, one of the candidates
for the medical code S46.211D in ICD10-CM is S46.212A, because these codes
share a common parent. Figure 1 presents examples of questions for the same
ICD10-CM code across each difficulty level. Algorithm 1 presents the pseudo-
code we used to generate the questions based on the required distance for each
level. We used this algorithm to create the medium and hard questions. Easy
questions were generated by a uniform sampling of the codes across the entire
codes’ graph.

Algorithm 1 Generating Questions by Distance

1:function
createQuestionsByDistance($\text{vocab\\_graph},\text{required\\_edges\\_distance},\text{num\\_of\\_options}$)

2: Input: Vocabulary graph vocab_graph, Required edges distance:
required_edges_distance, Number of options: num_of_options

3: Output: List of generated questions

4: $\text{nodes}\leftarrow\text{vocab\\_graph.nodes}$

5: $\text{questions}\leftarrow\varnothing$

6: for node in nodes do

7: $\text{candidate\\_nodes}\leftarrow\varnothing$

8: for candidate_node in $\text{nodes}\setminus\\{\text{node}\\}$ do

9:
$\text{path\\_length}\leftarrow\text{get\\_shortest\\_path\\_length}(\text{node},\text{candidate\\_node})$

10: if path_length in required_edges_distance then

11:
$\text{candidate\\_nodes}\leftarrow\text{candidate\\_nodes}\cup\\{\text{candidate\\_node}\\}$

12: end if

13: end for

14: $\text{candidates\\_number}\leftarrow\text{length of candidate\\_nodes}$

15: if $\text{candidates\\_number}\geq\text{num\\_of\\_options}$ then

16: $\text{options}\leftarrow\text{choose\\_random(candidate\\_nodes,
num\\_of\\_options)}$

17: $\text{question}\leftarrow\text{generate\\_question(node, options)}$

18: $\text{questions}\leftarrow\text{questions}\cup\\{\text{question}\\}$

19: end if

20: end for

21: return questions

22:end function

Figure 1: Questions across three difficulty levels for the S46.211D medical
code within the ICD10-CM vocabulary.

##  3 Experiments

In this section, we evaluate CLLMs on the MedConceptsQA benchmark. The CLLMs
we use for evaluations are: BioMistral (BioMistral-7B-DARE) [4], BioMedGPT
(BioMedGPT-LM-7B) [5], Gatortron (gatortron-large) [6], Llama3-OpenBioLLM
(70B) [20], BioBert [21], Meerket (7B) [7], Meditron (70B) [3], Clinical-
Longformer [22]. Additionally, we include GPT-3.5 and GPT-4 [17] for our
evaluation, although the main focus of these models is to be general-purpose
LLMs and not specifically focused on a clinical domain.

Our evaluation is based on zero-shot learning and few-shot learning for all
the models. For zero-shot learning we contain the question and general
instruction description. For few-shot learning we include 4 shots in the
prompt and then the question we want to get the answer for. Figure 2(a) shows
an example for zero-shot learning prompt for our drugs dataset (ATC) in the
benchmark, and figure 2(b) shows an example for few-shot learning example from
our benchmark, for ICD10-CM vocabulary.

###  3.1 Experiments Setup

For each dataset in the benchmark, we conduct zero-shot and few-shot
evaluations at each difficulty level (easy, medium, and hard). We repeat each
evaluation for each model three times and calculate a 95% confidence interval.
Due to the large amount of resources required by some of the models,
especially GPT4, we limit each type of test (vocabulary and difficulty and
shots) to 250 randomly sampled Q&As. Few-shot learning is performed using
4-shot learning. All models except GPT-3.5 (gpt-3.5-turbo) and GPT-4
(gpt-4-0125-preview) were evaluated using HuggingFace [23] on RTX-6000 GPU,
while GPT-3.5 and GPT-4 were inferred using the OpenAI API because these
models are not provided as open-source. We use accuracy as the evaluation
metric because the datasets in the benchmark are balanced, as we selected the
placement of the correct answer randomly during the creation of the benchmark.

(a) Zero-shot learning prompt example, for ATC vocabulary.

(b) Few-shot learning prompt example, for ICD10-CM vocabulary. We used 4
shots, each one contains the question and the answer.

Figure 2: Examples of zero-shot and few-shot learning prompts.

###  3.2 Aggregated Analysis

We perform zero-shot and few-shot evaluations for each model across all
vocabularies and difficulty levels (easy, medium, hard). In this section, we
present aggregated results for each model. The aggregation includes averaging
both the accuracy and the 95% confidence interval. Figure 3 presents the
results of zero-shot and few-shot learning evaluation. The full results across
all the models, vocabularies, and difficulty levels can be found in the
appendix (Table 3 for zero-shot learning and Table 3 for few-shot learning).

The CLLMs achieve performance close to Random guessing. For zero-shot learning
evaluation, None of the CLLMs outperforms the Random guessing threshold when
considering the error rate (confidence intervals). However, both GPT-3.5 and
GPT-4 outperform the CLLMs, with GPT-4 achieving the highest performance at
52.49% accuracy. This result represents a significant improvement over random
guessing, with an absolute increase of 27.49%. In the case of few-shot
learning, the results are similar for the CLLMs. However, GPT-3.5 shows an
improvement in accuracy of 4.418% (absolute), and GPT-4 demonstrates absolute
accuracy improvements of 9.422% between zero-shot and few-shot evaluations.

(a) Zero-shot aggregated results for each model.

(b) Few-shot aggregated results for each model.

Figure 3: Aggregated results for zero-shot and few-shot evaluations. The
vertical line represents the accuracy of random guessing for comparison.

###  3.3 Results for Difficulty Level

(a) Easy level, zero-shot learning.

(b) Easy level, few-shot learning.

(c) Medium level, zero-shot learning.

(d) Medium level, few-shot learning.

(e) Hard level, zero-shot learning.

(f) Hard level, few-shot learning.

Figure 4: Zero-shot and few-shot results for each of the levels (easy, medium,
and hard) with 95% confidence intervals over three runs. Results are
aggregated over difficulty vocabularies.

We create three difficulty levels (easy, medium, and hard) for each dataset in
our benchmark. Figure 4 shows the results of all models across the different
levels for both zero-shot and few-shot learning settings. Nearly all the CLLMs
achieve random guessing accuracy across levels (considering error bars),
except for LLaMA-3B OpenBioLLM-70B, which performs better than random guessing
on the easy level under few-shot learning. GPT-3.5 and GPT-4 demonstrate
improvement from zero-shot to few-shot learning, indicating that the 4-shot
examples enhance their performance. Moreover, the performance of GPT-3.5 and
GPT-4 deteriorates as task difficulty increases, validating that our benchmark
levels indeed represent a gradual increase in complexity. In Section 3.5, we
present a case study focusing on the ICD9-CM vocabulary and analyze results
across difficulty levels.

###  3.4 Results for Different Vocabularies

The accuracy results for each vocabulary (ATC, ICD9-CM, ICD10-CM, ICD9-PROC,
ICD10-PROC) are presented in Figure 5. The aggregation for the accuracy and
95% CI is for all the levels for each vocabulary. Also here, GPT-4
consistently outperform all other models across all vocabularies, with GPT-3.5
in second place. The CLLMs achieve performance close to random guessing. For
ATC, GPT-4 get an accuracy of 45.6% with zero-shot learning but 62.045% with
few-shot learning. In the case of ICD10-CM, GPT-4 demonstrates the highest
accuracy at 73.022% with few-shot learning, while for ICD9-CM, it achieves
74.711%. For procedures medical codes (ICD9-PROC and ICD10-PROC), the
performance is lower than the other vocabularies. The best model is also GPT-4
for both few-shot evaluations, achieving scores of 44.489% for ICD10-PROC and
55.289% for ICD9-PROC. Overall, these experiments show that almost all CLLMs
achieve performance close to random guessing, not only in the overall
aggregated score but also for each individual vocabulary. However, for ICD9-CM
and ICD10-CM, Llama3-OpenBioLLM-70B outperformed random guessing.

(a) ATC vocabulary, zero-shot

(b) ATC vocabulary, few-shot

(c) ICD10-CM vocabulary, zero-shot

(d) ICD10-CM vocabulary, few-shot

(e) ICD9-CM vocabulary, zero-shot

(f) ICD9-CM vocabulary, few-shot

(g) ICD10-PROC vocabulary, zero-shot

(h) ICD10-PROC vocabulary, few-shot

(i) ICD9-PROC vocabulary, zero-shot

(j) ICD9-PROC vocabulary, few-shot

Figure 5: Zero-shot and few-shot results for each of the vocabularies with 95%
confidence intervals over three runs. Results are aggregated over difficulty
levels.

###  3.5 Case Study: ICD9-CM Performance by Difficulty Level

Table 2 illustrates a case study of the ICD9-CM vocabulary and the performance
of the models across these levels (easy, medium, and hard). The objective of
these experiments is to highlight the diverse range of difficulty levels
represented by the questions within our proposed benchmark. As depicted in
this table, the results of both GPT-4 and GPT-3.5 decrease as the level
becomes more challenging. For easy questions, GPT-4 achieves very high
accuracy scores, with 98% for few-shot learning and 95.333% for zero-shot
learning. However, for medium questions, which are more difficult than the
easy ones, GPT-4 achieves only 65.467% for zero-shot and 72.4% for few-shot
learning. As for the hard level questions, the accuracy decreases
dramatically, with only 47.467% for zero-shot and 53.733% for few-shot
learning. This represents an absolute decrease of 44.267% in few-shot accuracy
between the easy and hard questions, and 47.866% for zero-shot.

Table 2: Zero-shot and Few-shot learning results with 95% confidence
intervals, for ICD9-CM vocabulary.

Model | Zero-shot | Few-shot
---|---|---
Easy | Medium | Hard | Easy | Medium | Hard
BioMistral | $25.333\pm 2.144$ | $24.933\pm 4.091$ | $27.333\pm 2.045$ | $25.6\pm 2.759$ | $23.067\pm 1.593$ | $26.8\pm 2.357$
BioBert | $23.333\pm 2.237$ | $27.467\pm 5.257$ | $26.133\pm 3.343$ | $25.467\pm 6.81$ | $25.067\pm 1.833$ | $25.067\pm 1.458$
Llama3–OpenBioLLM - 70B | $31.2\pm 7.298$ | $27.467\pm 9.174$ | $23.467\pm 4.114$ | $36.4\pm 0.454$ | $26.533\pm 3.666$ | $27.2\pm 1.571$
BioMedGPT | $24.0\pm 2.759$ | $25.067\pm 2.579$ | $22.133\pm 0.944$ | $24.533\pm 5.412$ | $24.267\pm 3.435$ | $24.267\pm 2.046$
Clinical-Longformer | $25.2\pm 0.454$ | $24.933\pm 3.404$ | $24.667\pm 2.328$ | $28.4\pm 1.636$ | $26.533\pm 4.808$ | $26.4\pm 1.2$
GPT-3.5 | $64.0\pm 0.454$ | $36.133\pm 1.309$ | $40.533\pm 1.593$ | $77.067\pm 3.282$ | $40.267\pm 3.282$ | $39.467\pm 1.718$
GPT-4 | $\textbf{95.333}\pm 2.328$ |  $\textbf{65.467}\pm$2.283 | $\textbf{47.467}\pm 4.214$ | $\textbf{98.0}\pm 1.2$ | $\textbf{72.4}\pm 1.815$ | $\textbf{53.733}\pm 2.499$
Gatortron | $24.667\pm 4.763$ | $26.933\pm 4.806$ | $26.133\pm 2.498$ | $25.6\pm 2.078$ | $26.0\pm 3.542$ | $24.667\pm 3.464$
MedMNX-7B | $23.6\pm 5.499$ | $24.4\pm 3.543$ | $22.8\pm 4.374$ | $24.267\pm 1.458$ | $25.067\pm 1.593$ | $25.2\pm 2.975$
Meditron-70B | $27.733\pm 4.655$ | $25.2\pm 2.4$ | $26.533\pm 3.694$ | $27.067\pm 5.316$ | $24.133\pm 4.634$ | $22.0\pm 5.684$
Meerkat | $25.2\pm 2.079$ | $27.467\pm 3.858$ | $25.333\pm 2.498$ | $26.667\pm 2.499$ | $27.067\pm 0.693$ | $24.133\pm 1.458$

##  4 Discussion

(a) Meditron-70B error analysis.

(b) Meerkat-7B error analysis.

(c) BioMistral-7B error analysis.

Figure 6: Error analysis (confusion matrix) of few-shot learning for ICD9-CM
at the easy level, for three CLLMs that achieved random guessing performance
on our benchmark.

With the proposed benchmark in this work, we evaluate the understanding of
medical concepts across LLMs, including those not pre-trained on medical data.
Our experiments show that the current state-of-the-art CLLMs are not adequate
in understanding the meaning of medical concepts such as diagnoses,
procedures, and drug codes, and the differences between them. The performances
of these CLLMs were close to random guessing. However, GPT-3.5 and GPT-4
achieved better accuracy for our benchmark, although the accuracy is still not
high enough in some cases. For example, the few-shot learning accuracy of
GPT-4 for ICD10-PROC and ICD9-PROC is only 44.489% and 55.289%, respectively.
Additionally, we present a case study of the ICD9-CM vocabulary, which shows a
significant performance difference between GPT-4 and GPT-3.5 models across
different difficulty levels, demonstrating that our benchmark indeed contains
questions of varying difficulty.

In order to investigate why the CLLMs achieve a random guessing guess
performance, we calculate the confusion matrix for three CLLMs (Meditron-70B,
Meerkat-7B, and BioMistral-7B) for our case study vocabulary, ICD9-CM, under
few-shot learning settings. We choose this level because it is the easiest
level on the benchmark, and GPT-4 achieved very high accuracy, however, the
CLLMs also achieved Random guessing performance on this level. Figure 6 shows
the error analysis using three confusion matrices for Meditron-70B,
BioMistral-7B-DARE, and Meerkat-7B on the ICD9-CM easy level. Meditron-70B
almost always predicts the label B, and Meerkat-7B almost always predicts A or
D but never predicts C. BioMistral-7B does not always predict the same label
but predicts diverse answers among the given options (A, B, C, D). However,
the answers from this model are mostly incorrect since it also achieved Random
guessing accuracy. Our careful dataset design, incorporating diverse label
representation through random sampling in the 4-shot examples, ensures that we
did not inadvertently influence the models (Meditron-70B and Meerkat-7B) to
consistently predict the same label for this dataset. We encourage further
studies to investigate this issue, focusing on how effectively existing CLLMs
follow instructions, whether in zero-shot or few-shot learning scenarios.

##  5 Limitations and Future Work

Limitations: We did not include all medical code vocabularies, for example,
LOINC, which is a vocabulary for identifying medical laboratory observations.
Additionally, our benchmark contains questions in only a single form.

Future Work: In this work, we demonstrated that CLLMs achieved random guessing
performance on the MedConceptsQA benchmark. Future research should focus on
improving CLLMs to attain a better understanding of medical concepts, aiming
for performance comparable to or better than GPT-4, whose primary focus is not
medical, unlike CLLMs. Additionally, since our benchmark is competitive and
evaluates a LLMs understanding of medical concepts, we encourage new CLLMs to
include MedConceptsQA in their evaluation, with the aim of achieving superior
results compared to current CLLMs.

##  6 Conclusion

In this study, we introduced MedConceptsQA, an open-source benchmark
comprising over 800,000 questions spanning three difficulty levels. It was
designed to evaluate the understanding and reasoning capabilities of LLMs on
medical concepts across diagnoses, procedures, and drugs. Experimental results
showed that state-of-the-art CLLMs, despite being pre-trained on medical data,
achieved accuracy levels close to random guessing on this benchmark. However,
general-purpose models (GPT-3.5 and GPT-4) outperformed CLLMs. Notably, GPT-4
exhibited the best performance, although its accuracy remained insufficient
for certain datasets in our benchmark.

##  7 Reproducibility

The MedConceptsQA benchmark is available at
https://huggingface.co/datasets/ofir408/MedConceptsQA.

Our code for creating this benchmark, and the evaluation code are available at
the following link: https://github.com/nadavlab/MedConceptsQA.

## References

  * [1] Wayne Xin Zhao et al.  “A survey of large language models”  In _arXiv preprint arXiv:2303.18223_ , 2023
  * [2] Jean Kaddour et al.  “Challenges and applications of large language models”  In _arXiv preprint arXiv:2307.10169_ , 2023
  * [3] Zeming Chen et al.  “Meditron-70b: Scaling medical pretraining for large language models”  In _arXiv preprint arXiv:2311.16079_ , 2023
  * [4] Yanis Labrak et al.  “BioMistral: A Collection of Open-Source Pretrained Large Language Models for Medical Domains”  In _arXiv preprint arXiv:2402.10373_ , 2024
  * [5] Yizhen Luo et al.  “Biomedgpt: Open multimodal generative pre-trained transformer for biomedicine”  In _arXiv preprint arXiv:2308.09442_ , 2023
  * [6] Xi Yang et al.  “A large language model for electronic health records”  In _npj Digital Medicine_ 5.1 Nature Publishing Group UK London, 2022, pp. 194
  * [7] Hyunjae Kim et al.  “Small Language Models Learn Enhanced Reasoning Skills from Medical Textbooks”  In _arXiv preprint arXiv:2404.00376_ , 2024
  * [8] Ofir Ben Shoham and Nadav Rappoport  “Cpllm: Clinical prediction with large language models”  In _arXiv preprint arXiv:2309.11295_ , 2023
  * [9] Yuting He et al.  “Foundation Model for Advancing Healthcare: Challenges, Opportunities, and Future Directions”  In _arXiv preprint arXiv:2404.03264_ , 2024
  * [10] Anastasia Krithara, Anastasios Nentidis, Konstantinos Bougiatiotis and Georgios Paliouras  “BioASQ-QA: A manually curated corpus for Biomedical Question Answering”  In _Scientific Data_ 10.1 Nature Publishing Group UK London, 2023, pp. 170
  * [11] Anusri Pampari, Preethi Raghavan, Jennifer Liang and Jian Peng  “emrQA: A Large Corpus for Question Answering on Electronic Medical Records”  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , 2018, pp. 2357–2368
  * [12] Qiao Jin et al.  “PubMedQA: A Dataset for Biomedical Research Question Answering”  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ Hong Kong, China: Association for Computational Linguistics, 2019, pp. 2567–2577  DOI: 10.18653/v1/D19-1259
  * [13] Olumide E Ojo et al.  “MedAI Dialog Corpus (MEDIC): Zero-Shot Classification of Doctor and AI Responses in Health Consultations”  In _arXiv preprint arXiv:2310.12489_ , 2023
  * [14] Diego Mollá, María Elena Santiago-Martínez, Abeed Sarker and Cécile Paris  “A corpus for research in text processing for evidence based medicine”  In _Language Resources and Evaluation_ 50 Springer, 2016, pp. 705–727
  * [15] Amin Dada et al.  “CLUE: A Clinical Language Understanding Evaluation for LLMs”  In _arXiv preprint arXiv:2404.04067_ , 2024
  * [16] Ali Soroush et al.  “Large Language Models Are Poor Medical Coders—Benchmarking of Medical Code Querying”  In _NEJM AI_ Massachusetts Medical Society, 2024, pp. AIdbp2300040
  * [17] Josh Achiam et al.  “Gpt-4 technical report”  In _arXiv preprint arXiv:2303.08774_ , 2023
  * [18] Shengyu Zhang et al.  “Instruction tuning for large language models: A survey”  In _arXiv preprint arXiv:2308.10792_ , 2023
  * [19] Chaoqi Yang et al.  “PyHealth: A Deep Learning Toolkit for Healthcare Predictive Modeling”  In _Proceedings of the 27th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD) 2023_ , 2023  URL: https://github.com/sunlabuiuc/PyHealth
  * [20] Malaikannan Sankarasubbu Ankit Pal  “OpenBioLLMs: Advancing Open-Source Large Language Models for Healthcare and Life Sciences”  In _Hugging Face repository_ Hugging Face, https://huggingface.co/aaditya/OpenBioLLM-Llama3-70B, 2024
  * [21] Jinhyuk Lee et al.  “BioBERT: a pre-trained biomedical language representation model for biomedical text mining”  In _Bioinformatics_ 36.4 Oxford University Press, 2020, pp. 1234–1240
  * [22] Yikuan Li et al.  “A comparative study of pretrained language models for long clinical text”  In _Journal of the American Medical Informatics Association_ 30.2 Oxford University Press, 2023, pp. 340–347
  * [23] Thomas Wolf et al.  “Huggingface’s transformers: State-of-the-art natural language processing”  In _arXiv preprint arXiv:1910.03771_ , 2019

##  8 Appendix

Table 3 presents the comprehensive zero-shot learning results for all models,
vocabularies, and difficulty levels (easy, medium, and hard). Table 4 presents
the few-shot learning results.

| ICD9PROC | ICD9CM | ICD10PROC | ICD10CM | ATC
---|---|---|---|---|---
| easy | medium | hard | easy | medium | hard | easy | medium | hard | easy | medium | hard | easy | medium | hard
model |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
BioMistral/BioMistral-7B-DARE | $22.4\pm 2.08$ | $24.8\pm 2.4$ | $25.2\pm 3.54$ | $25.33\pm 2.14$ | $24.93\pm 4.09$ | $27.33\pm 2.05$ | $22.13\pm 3.52$ | $19.47\pm 5.22$ | $26.0\pm 4.08$ | $23.07\pm 4.54$ | $24.27\pm 2.88$ | $25.47\pm 2.58$ | $22.8\pm 2.97$ | $31.33\pm 12.31$ | $24.0\pm 3.6$
HuggingFaceH4/zephyr-7b-beta | $23.6\pm 2.4$ | $23.2\pm 3.18$ | $26.67\pm 0.52$ | $25.87\pm 4.54$ | $27.07\pm 2.33$ | $26.0\pm 1.64$ | $25.2\pm 3.54$ | $25.33\pm 3.02$ | $26.0\pm 6.68$ | $26.0\pm 2.27$ | $26.13\pm 4.36$ | $24.27\pm 1.05$ | $26.4\pm 3.63$ | $24.53\pm 3.02$ | $26.8\pm 3.95$
PharMolix/BioMedGPT-LM-7B | $24.0\pm 3.63$ | $23.73\pm 6.02$ | $23.07\pm 2.14$ | $24.0\pm 2.76$ | $25.07\pm 2.58$ | $22.13\pm 0.94$ | $24.53\pm 5.22$ | $29.47\pm 9.70$ | $27.47\pm 2.28$ | $22.8\pm 2.72$ | $26.53\pm 4.47$ | $25.6\pm 2.4$ | $25.87\pm 3.43$ | $21.73\pm 9.99$ | $25.2\pm 4.99$
UFNLP/gatortron-large | $25.47\pm 3.67$ | $27.73\pm 1.39$ | $22.8\pm 0.45$ | $24.67\pm 4.76$ | $26.93\pm 4.81$ | $26.13\pm 2.50$ | $21.87\pm 5.63$ | $24.53\pm 2.33$ | $25.33\pm 1.89$ | $25.87\pm 2.33$ | $23.87\pm 2.77$ | $26.53\pm 3.28$ | $23.6\pm 4.08$ | $19.6\pm 3.42$ | $23.07\pm 3.18$
aaditya/Llama3-OpenBioLLM-70B | $25.87\pm 4.89$ | $28.13\pm 9.30$ | $24.53\pm 4.12$ | $31.2\pm 7.30$ | $27.47\pm 9.18$ | $23.47\pm 4.12$ | $25.47\pm 3.18$ | $27.6\pm 7.20$ | $22.8\pm 7.87$ | $29.73\pm 3.86$ | $26.27\pm 4.54$ | $25.73\pm 0.69$ | $27.33\pm 4.59$ | $26.0\pm 5.52$ | $22.4\pm 1.2$
dmis-lab/biobert-v1.1 | $25.2\pm 4.54$ | $25.2\pm 3.18$ | $27.73\pm 3.43$ | $23.33\pm 2.24$ | $27.47\pm 5.26$ | $26.13\pm 3.34$ | $26.27\pm 2.50$ | $24.0\pm 5.52$ | $24.93\pm 4.54$ | $27.6\pm 1.64$ | $28.13\pm 5.06$ | $25.33\pm 2.58$ | $26.27\pm 4.12$ | $26.0\pm 2.72$ | $28.67\pm 2.92$
dmis-lab/meerkat-7b-v1.0 | $24.67\pm 3.67$ | $25.2\pm 6.24$ | $22.4\pm 4.74$ | $25.2\pm 2.08$ | $27.47\pm 3.86$ | $25.33\pm 2.50$ | $23.73\pm 1.39$ | $29.73\pm 15.98$ | $27.6\pm 2.53$ | $24.67\pm 1.14$ | $28.67\pm 1.39$ | $25.47\pm 6.69$ | $25.2\pm 2.97$ | $24.0\pm 2.76$ | $25.87\pm 1.83$
epfl-llm/meditron-70b | $24.53\pm 6.04$ | $26.67\pm 5.06$ | $24.4\pm 1.36$ | $27.73\pm 4.65$ | $25.2\pm 2.4$ | $26.53\pm 3.69$ | $24.93\pm 1.31$ | $27.87\pm 2.10$ | $22.53\pm 0.52$ | $25.87\pm 1.05$ | $23.07\pm 0.52$ | $25.6\pm 0.79$ | $27.87\pm 0.52$ | $23.73\pm 2.62$ | $23.87\pm 6.81$
epfl-llm/meditron-7b | $25.33\pm 3.86$ | $25.87\pm 2.05$ | $27.2\pm 2.83$ | $26.53\pm 0.69$ | $24.93\pm 2.28$ | $25.07\pm 2.50$ | $23.07\pm 4.98$ | $25.47\pm 4.07$ | $24.8\pm 1.36$ | $25.47\pm 1.72$ | $24.67\pm 9.13$ | $28.27\pm 3.80$ | $26.0\pm 2.4$ | $27.73\pm 5.16$ | $25.87\pm 3.28$
gpt-3.5-turbo | $31.2\pm 2.76$ | $33.33\pm 1.14$ | $29.07\pm 3.67$ | $64.0\pm 0.45$ | $36.13\pm 1.31$ | $40.53\pm 1.59$ | $24.0\pm 2.97$ | $35.87\pm 3.67$ | $28.53\pm 2.05$ | $68.27\pm 3.02$ | $37.6\pm 4.80$ | $37.07\pm 2.05$ | $35.33\pm 2.73$ | $25.33\pm 0.94$ | $29.6\pm 2.83$
gpt-4-0125-preview | $54.0\pm 6.28$ | $40.13\pm 4.09$ | $35.73\pm 0.94$ | $95.33\pm 2.33$ | $65.47\pm 2.28$ | $47.47\pm 4.21$ | $30.0\pm 3.14$ | $42.27\pm 3.43$ | $34.93\pm 1.14$ | $94.27\pm 2.14$ | $54.53\pm 4.54$ | $56.4\pm 2.4$ | $57.87\pm 3.02$ | $45.33\pm 4.98$ | $33.6\pm 2.08$
johnsnowlabs/JSL-MedMNX-7B | $22.67\pm 3.22$ | $23.47\pm 1.39$ | $22.53\pm 2.58$ | $23.6\pm 5.50$ | $24.4\pm 3.54$ | $22.8\pm 4.37$ | $27.07\pm 1.72$ | $23.6\pm 3.54$ | $26.27\pm 2.10$ | $25.47\pm 3.02$ | $24.67\pm 4.07$ | $24.4\pm 1.2$ | $26.13\pm 1.46$ | $24.93\pm 4.29$ | $24.4\pm 5.79$
meta-llama/Meta-Llama-3-8B-Instruct | $24.0\pm 1.57$ | $28.27\pm 6.93$ | $27.73\pm 5.83$ | $26.67\pm 8.91$ | $25.87\pm 3.46$ | $26.67\pm 3.64$ | $22.13\pm 2.92$ | $32.67\pm 5.22$ | $22.0\pm 9.47$ | $27.6\pm 13.75$ | $22.67\pm 6.42$ | $24.8\pm 9.81$ | $26.67\pm 5.67$ | $25.2\pm 5.52$ | $24.67\pm 3.86$
yikuan8/Clinical-Longformer | $26.4\pm 2.4$ | $27.2\pm 1.2$ | $24.0\pm 5.11$ | $25.2\pm 0.45$ | $24.93\pm 3.40$ | $24.67\pm 2.33$ | $24.4\pm 3.6$ | $22.93\pm 3.78$ | $24.27\pm 2.33$ | $28.4\pm 2.53$ | $24.0\pm 2.08$ | $24.0\pm 1.36$ | $25.73\pm 1.59$ | $23.2\pm 0.91$ | $26.27\pm 3.02$

Table 3: Zero-shot learning full results for all the models, vocabularies, and
levels.

| ICD9PROC | ICD9CM | ICD10PROC | ICD10CM | ATC
---|---|---|---|---|---
| easy | medium | hard | easy | medium | hard | easy | medium | hard | easy | medium | hard | easy | medium | hard
model |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
BioMistral/BioMistral-7B-DARE | 29.467 $\pm$ 2.045 | 23.867 $\pm$ 3.694 | 25.6 $\pm$ 2.759 | 25.6 $\pm$ 2.759 | 23.067 $\pm$ 1.593 | 26.8 $\pm$ 2.356 | 26.8 $\pm$ 2.078 | 22.0 $\pm$ 2.078 | 21.067 $\pm$ 3.91 | 24.0 $\pm$ 3.6 | 26.0 $\pm$ 2.356 | 22.8 $\pm$ 3.142 | 24.8 $\pm$ 0.907 | 27.867 $\pm$ 5.707 | 26.133 $\pm$ 1.141
HuggingFaceH4/zephyr-7b-beta | 25.867 $\pm$ 3.02 | 23.733 $\pm$ 3.434 | 26.8 $\pm$ 0.454 | 24.133 $\pm$ 1.047 | 26.933 $\pm$ 6.257 | 25.6 $\pm$ 1.977 | 24.8 $\pm$ 2.759 | 23.2 $\pm$ 0.785 | 25.333 $\pm$ 2.045 | 23.333 $\pm$ 0.944 | 27.333 $\pm$ 2.498 | 26.533 $\pm$ 2.143 | 23.6 $\pm$ 1.2 | 24.933 $\pm$ 2.095 | 23.733 $\pm$ 1.141
PharMolix/BioMedGPT-LM-7B | 23.467 $\pm$ 2.916 | 29.733 $\pm$ 7.539 | 24.667 $\pm$ 1.833 | 24.533 $\pm$ 5.41 | 24.267 $\pm$ 3.434 | 24.267 $\pm$ 2.045 | 24.4 $\pm$ 1.2 | 22.533 $\pm$ 3.962 | 24.667 $\pm$ 5.67 | 24.533 $\pm$ 1.593 | 22.4 $\pm$ 2.4 | 23.6 $\pm$ 1.635 | 26.267 $\pm$ 3.281 | 26.133 $\pm$ 4.543 | 28.4 $\pm$ 2.974
UFNLP/gatortron-large | 25.2 $\pm$ 3.6 | 23.2 $\pm$ 2.525 | 24.533 $\pm$ 4.358 | 25.6 $\pm$ 2.078 | 26.0 $\pm$ 3.542 | 24.667 $\pm$ 3.464 | 25.067 $\pm$ 1.888 | 25.2 $\pm$ 3.175 | 23.867 $\pm$ 2.951 | 24.8 $\pm$ 1.635 | 24.0 $\pm$ 3.542 | 24.933 $\pm$ 0.693 | 26.667 $\pm$ 3.217 | 23.333 $\pm$ 2.283 | 25.867 $\pm$ 0.944
aaditya/Llama3-OpenBioLLM-70B | 29.467 $\pm$ 3.217 | 25.333 $\pm$ 0.693 | 22.533 $\pm$ 1.141 | 36.4 $\pm$ 0.454 | 26.533 $\pm$ 3.666 | 27.2 $\pm$ 1.571 | 24.667 $\pm$ 2.951 | 26.267 $\pm$ 1.309 | 28.667 $\pm$ 2.951 | 31.867 $\pm$ 3.02 | 26.8 $\pm$ 3.6 | 24.8 $\pm$ 1.977 | 27.067 $\pm$ 3.343 | 23.067 $\pm$ 3.694 | 22.4 $\pm$ 0.785
dmis-lab/biobert-v1.1 | 28.533 $\pm$ 3.343 | 23.867 $\pm$ 1.141 | 23.867 $\pm$ 2.327 | 25.467 $\pm$ 6.808 | 25.067 $\pm$ 1.833 | 25.067 $\pm$ 1.458 | 26.667 $\pm$ 1.385 | 25.6 $\pm$ 3.424 | 26.267 $\pm$ 1.047 | 24.8 $\pm$ 2.4 | 27.333 $\pm$ 3.087 | 27.333 $\pm$ 3.91 | 26.0 $\pm$ 4.799 | 22.4 $\pm$ 1.977 | 23.6 $\pm$ 0.785
dmis-lab/meerkat-7b-v1.0 | 26.4 $\pm$ 5.948 | 24.533 $\pm$ 1.717 | 24.933 $\pm$ 3.694 | 26.667 $\pm$ 2.498 | 27.067 $\pm$ 0.693 | 24.133 $\pm$ 1.458 | 24.267 $\pm$ 1.141 | 27.6 $\pm$ 7.71 | 23.867 $\pm$ 2.579 | 26.933 $\pm$ 2.498 | 22.133 $\pm$ 3.666 | 23.733 $\pm$ 4.065 | 23.467 $\pm$ 2.045 | 23.467 $\pm$ 2.143 | 24.933 $\pm$ 1.309
epfl-llm/meditron-70b | 25.2 $\pm$ 2.078 | 25.867 $\pm$ 0.944 | 26.533 $\pm$ 3.281 | 27.067 $\pm$ 5.315 | 24.133 $\pm$ 4.632 | 22.0 $\pm$ 5.682 | 27.067 $\pm$ 2.045 | 24.8 $\pm$ 3.927 | 25.333 $\pm$ 1.888 | 25.867 $\pm$ 3.962 | 26.133 $\pm$ 7.045 | 25.867 $\pm$ 4.72 | 24.533 $\pm$ 2.143 | 24.4 $\pm$ 2.974 | 24.133 $\pm$ 1.833
epfl-llm/meditron-7b | 24.4 $\pm$ 1.361 | 26.133 $\pm$ 3.803 | 25.6 $\pm$ 0.907 | 24.4 $\pm$ 2.721 | 22.8 $\pm$ 4.031 | 24.0 $\pm$ 1.2 | 25.467 $\pm$ 2.771 | 22.267 $\pm$ 9.991 | 23.733 $\pm$ 4.214 | 21.867 $\pm$ 4.632 | 24.133 $\pm$ 3.637 | 22.533 $\pm$ 3.02 | 24.667 $\pm$ 2.657 | 21.067 $\pm$ 5.157 | 23.733 $\pm$ 2.327
gpt-3.5-turbo | 31.733 $\pm$ 6.09 | 31.067 $\pm$ 2.045 | 30.8 $\pm$ 0.907 | 77.067 $\pm$ 3.281 | 40.267 $\pm$ 3.281 | 39.467 $\pm$ 1.717 | 27.333 $\pm$ 0.262 | 28.533 $\pm$ 2.771 | 29.067 $\pm$ 2.771 | 76.0 $\pm$ 4.326 | 39.6 $\pm$ 1.635 | 39.067 $\pm$ 2.657 | 67.067 $\pm$ 0.944 | 39.2 $\pm$ 1.571 | 25.867 $\pm$ 2.951
gpt-4-0125-preview | 80.0 $\pm$ 4.467 | 50.8 $\pm$ 1.814 | 35.067 $\pm$ 3.434 | 98.0 $\pm$ 1.2 | 72.4 $\pm$ 1.814 | 53.733 $\pm$ 2.498 | 46.8 $\pm$ 3.424 | 52.267 $\pm$ 1.458 | 34.4 $\pm$ 1.571 | 98.133 $\pm$ 0.944 | 63.867 $\pm$ 2.498 | 57.067 $\pm$ 2.88 | 92.4 $\pm$ 0.785 | 59.467 $\pm$ 4.543 | 34.267 $\pm$ 1.458
johnsnowlabs/JSL-MedMNX-7B | 25.067 $\pm$ 0.262 | 25.6 $\pm$ 2.356 | 28.133 $\pm$ 3.523 | 24.267 $\pm$ 1.458 | 25.067 $\pm$ 1.593 | 25.2 $\pm$ 2.974 | 24.133 $\pm$ 3.857 | 27.6 $\pm$ 1.2 | 24.667 $\pm$ 3.185 | 26.0 $\pm$ 3.175 | 27.333 $\pm$ 4.065 | 25.6 $\pm$ 2.356 | 22.4 $\pm$ 2.078 | 27.867 $\pm$ 4.214 | 25.467 $\pm$ 1.141
meta-llama/Meta-Llama-3-8B-Instruct | 27.6 $\pm$ 3.175 | 24.933 $\pm$ 5.41 | 25.6 $\pm$ 0.785 | 27.067 $\pm$ 1.833 | 24.267 $\pm$ 1.593 | 26.133 $\pm$ 1.141 | 24.8 $\pm$ 0.907 | 28.933 $\pm$ 4.995 | 26.133 $\pm$ 2.237 | 23.333 $\pm$ 0.944 | 21.867 $\pm$ 0.944 | 26.267 $\pm$ 4.654 | 26.133 $\pm$ 1.458 | 26.0 $\pm$ 8.276 | 25.733 $\pm$ 2.237
yikuan8/Clinical-Longformer | 26.0 $\pm$ 4.031 | 24.8 $\pm$ 2.525 | 26.4 $\pm$ 1.977 | 28.4 $\pm$ 1.635 | 26.533 $\pm$ 4.807 | 26.4 $\pm$ 1.2 | 23.867 $\pm$ 3.857 | 24.933 $\pm$ 0.944 | 25.2 $\pm$ 2.078 | 24.533 $\pm$ 2.498 | 26.4 $\pm$ 7.71 | 25.733 $\pm$ 5.448 | 24.267 $\pm$ 2.095 | 24.4 $\pm$ 5.517 | 25.333 $\pm$ 6.09

Table 4: Few-shot learning full results for all the models, vocabularies, and
levels.