# Explosions in Roche-Lobe Distorted Stars: Relativistic Bullets in Binaries

Marcus DuPont Center for Cosmology and Particle Physics, Physics Department,
New York University New York, NY, 10003, USA  Max-Planck-Institut für
Astrophysik, Karl-Schwarzschild-Stra$\beta$e 1, D-85741 Garching, Germany
Andrew MacFadyen Center for Cosmology and Particle Physics, Physics
Department, New York University New York, NY, 10003, USA  Selma de Mink Max-
Planck-Institut für Astrophysik, Karl-Schwarzschild-Stra$\beta$e 1, D-85741
Garching, Germany  Astronomical Institute Anton Pannekoek, Amsterdam
University, Science Park 904, 1098 XH, Amsterdam, Netherlands

###### Abstract

State-of-the-art surveys reveal that most massive stars in the universe evolve
in close binaries. Massive stars in such systems are expected to develop
aspherical envelopes due to tidal interactions and/or rotational effects.
Recently, it was shown that point explosions in oblate stars can produce
relativistic equatorial ring-like outflows. Moreover, since stripped-envelope
stars in binaries can expand enough to fill their Roche lobes anew, it is
likely that these stars die with a greater degree of asphericity than the
oblate spheroid geometry previously studied. We investigate the effects of
this asymmetry by studying the gas dynamics of axisymmetric point explosions
in stars in various stages of filling their Roche lobes. We find that point
explosions in these pear-shaped stars produce trans-relativistic ejecta that
coalesces into bullets pointed both toward and away from the binary companion.
We present this result and comment on key morphological differences between
core-collapse explosions in spherical versus distorted stars in binary
systems, effects on gravitational wave sources, and observational signatures
that could be used to glean these explosion geometries from current and future
surveys.

Roche lobe (1405) — Relativistic fluid dynamics (1389)— Binary stars (154)—
Supernovae (1668)— Radio Transient Sources (2008)

††software: CMasher (van der Velden, 2020), SIMBI (DuPont, 2023), MESA (Paxton
et al., 2011, 2013, 2015, 2018, 2019), astropy (Astropy Collaboration et al.,
2013)

##  1 Introduction

Massive stars predominantly form in binary and multiple systems and the
majority of them will interact and exchange mass with at least one companion
star (Sana et al., 2012; Moe & Di Stefano, 2017). Because these stars are so
massive, they will die as core-collapse supernovae (CCSNe), explosions that
synthesize many of the heavy elements and drive outflows from galaxies (e.g.
Fowler & Hoyle, 1964; Woosley et al., 2002; Efstathiou, 2000).

The shape of a star in a binary system can deviate significantly from
spherical symmetry. This is particularly true when a star expands to fill its
Roche lobe. This forces the star to become “pear-shaped” due to the combined
effect of the gravitational pull of the companion and the centrifugal effect
from the motion around the common center of mass. In this paper, we explore
what happens if a supernova explosion occurs in such a highly asymmetric,
Roche-lobe-deformed star.

Explosions in non-spherical stars are of interest because of the effects that
an oblique shock can have when breaking out of the stellar surface. The non-
radial motions can dramatically alter the dynamics and emission of shock
emergence (Matzner et al., 2013). In a previous paper, DuPont et al. (2022)
numerically explored the effect of point explosions in rapidly rotating,
oblate progenitors and showed how these produce relativistic equatorial ring-
like outflows. Here we build upon this work and extend this to the pear-shaped
stars in binary systems.

The fraction of core-collapse supernovae progenitors that are significantly
Roche-lobe deformed is not very well known, but binary evolutionary
simulations suggest that this may be common for lower mass progenitors of
stripped-envelope supernovae (type Ib/c and type IIb, Filippenko 1997). These
are produced in binary systems where the progenitor loses most of its H-rich
envelope during a first mass transfer phase (Kippenhahn & Weigert, 1967;
Götberg et al., 2017). The resulting stripped star first contracts as a hot
helium star. Drout et al. (2023) and Götberg et al. (2023) recently presented
the discovery and analysis of such stars, which so far had been missing.

During their later evolution, helium stars in the mass range $2-4M_{\odot}$
can expand significantly to become helium giants (e.g., Divine, 1965; Habets,
1986; Dewi et al., 2002; Dewi & Pols, 2003a; Yoon et al., 2010; Sravan et al.,
2019). If they retain a layer of hydrogen, the expansion is even more drastic
(Laplace et al., 2020; Klencki et al., 2022). This implies that these low-mass
stripped progenitors are (close to) filling their Roche lobes at the onset of
core collapse and therefore deformed.

Of special interest are supernova explosions in Roche-lobe deformed stripped
stars where the companion star is a neutron star. These systems may lead to
ultra-stripped SNe (Tauris et al., 2015). This is the last stage in the
pathway towards the formation of double neutron stars and gravitational wave
sources (Tauris et al., 2017). In this work, we present a suite of two-
dimensional axisymmetric simulations of point explosions in a $3M_{\odot}$
pre-supernova helium star at various stages of filling its Roche lobe,
described in Section 2. We study how the shock propagates and breaks out from
the surface. We report on the formation of a trans-relativistic ejecta pointed
both away and towards the companion in Section 3. We discuss the astrophysical
implications and observational consequences in Section 4. Section 5 provides a
summary with conclusions.

##  2 Numerical Setup

###  2.1 Governing Equations

The governing hydrodynamical equations in the special relativistic regime
describe the conservation of baryon number,

$\partial_{\mu}(\rho u^{\mu})=0,$ (1)

and conservation of stress-energy,

$\partial_{\nu}T^{\mu\nu}=0,$ (2)

where we work in the fluid frame so that $T^{\mu\nu}=\rho
hu^{\mu}u^{\nu}+p\eta^{\mu\nu}$. In these equations, $h=1+\epsilon+p/\rho$ is
the total specific enthalpy, $\epsilon$ is the specific internal energy,
$\rho$ is proper mass density, $\eta^{\mu\nu}$ is the Minkowski metric, $p$ is
pressure, $[u^{\mu}]=\Gamma(1,\vec{\beta})$ is four-velocity,
$\Gamma=(1-\beta^{2})^{-1/2}$ is Lorentz factor, and $\beta$ is the fluid
velocity in units of $c$. The signature of the metric is $(-+++)$. The system
of equations is closed by assuming an adiabatic equation of state
$p=(\hat{\gamma}-1)\rho\epsilon$, where $\hat{\gamma}=4/3$ is the adiabatic
gas index. We work in units where $c=1$. The point explosions are launched in
a 2D axisymmetric spherical-polar mesh in the $(r-\phi)$ plane. All quantities
are made dimensionless by combinations of $R_{\odot}$, $M_{\odot}$, and $c$.
The code used to run the simulations, entitled SIMBI, was written by this
study’s first author, and is an open-source, second-order, Godunov-type
scheme. SIMBI is is parallelized for multi-CPU frameworks enabled by OpenMP,
and supports GPU-acceleration on both AMD and Nvidia platforms (DuPont, 2023).

###  2.2 Initial Conditions

Using the MESA (Paxton et al., 2011, 2013, 2015, 2018, 2019) code, Farmer et
al. (2021) evolved a series of binary stars with a focus on the primary star
post helium burning. The authors followed the evolution of said primary stars
just before core collapse. The binary system had various initial masses, but
constant initial mass ratio $q_{0}=M_{2\rm init}/M_{\rm init}=0.8$, and each
star was non-rotating. We use these stars as the progenitors in this study.
From knowing the initial mass, $M_{\rm init}$, final mass, $M_{1}$, and
initial mass ratio, $q_{0}$, we can arrive at the final mass ratio of the
system pre-supernova assuming the mass exchange was conservative,

$q=M_{2}/M_{1}=M_{\rm init}(1+q_{0})/M_{1}-1.$ (3)

This parameter is crucial for constructing the equipotential surfaces of the
Roche lobe. The separation of the stars can be calculated by inverting the
Eggleton (1983) formula to arrive at

$d=R_{\rm eq}\left(\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln(1+q^{1/3})}\right)^{-1},$
(4)

where $R_{\rm eq}$ is the volume equivalent radius defined such that the
volume of the Roche lobe is $V=4/3\pi R_{\rm eq}^{3}$. For simplicity, we
choose a single primary star from the grid of binary simulations for our
entire analysis. This primary had an initial mass of $M_{\rm
init}=11M_{\odot}$ and pre-supernova mass of $M_{1}=3M_{\odot}$, which gives a
final mass ratio of $q=5.21$. Its volume equivalent radius is roughly
$5.4R_{\odot}$.

In this analysis, we assume a steady wind environment,

$\rho_{\rm amb}(r,\phi)=\frac{\dot{M}}{4\pi v_{\rm
wind}R_{*}^{2}(r,\phi)}\left(\frac{r}{R_{*}(r,\phi)}\right)^{-2}\equiv\frac{A}{r^{2}},$
(5)

where $\dot{M}$ is the stellar mass loss rate, $v_{\rm wind}$ is the velocity
of the wind, and $R_{*}$ is the outer radius of the distorted envelope. A
common convention for modulating stellar wind magnitudes is the mass-loading
parameter $A=A_{*}\times 5\times 10^{11}\,\mathrm{g\ cm^{-1}}$, which assumes
$\dot{M}=10^{-5}M_{\odot}\,\mathrm{yr^{-1}}$ and $v_{\rm
wind}=10^{8}\,\mathrm{cm\ s^{-1}}$. We fix $A_{*}=1$ in our simulations.

The MESA density profiles are computed in 1D using the volume equivalent
radius, $R_{\rm eq}$. Since the equipotential surfaces of the Roche lobe are
themselves isobars, we can invoke hydrostatic equilibrium ($-\nabla
P/\nabla\Phi=\rho$) to map the density $\rho(R_{\rm eq})$ to the surfaces of
constant potential, $\Phi$. We also assume synchronous rotation for the system
such that the dimensionless potential has the form (Kopal, 1959):

$\displaystyle\Phi(r,\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{r}+\frac{q+1}{2}r^{2}$ $\displaystyle+$ $\displaystyle
q\left(\frac{1}{\sqrt{1-2r\cos\phi+r^{2}}}-r\cos\phi\right),$

where we set $\theta=\pi/2$ so that $r-\phi$ lies in the orbital plane of the
binary. The formulae above allow us to transform the 1D density profiles from
$\rho(R_{\rm eq})$ to multidimensional profiles $\rho(\Phi)$, conserving the
stellar mass in the process. To map the density contours for a Roche lobe
filling factor, $\mathcal{F}$, we use the functional form of the modified
potential described in Leahy & Leahy (2015):

$\Phi_{\mathcal{F}}=\frac{\Phi_{\rm
L1}+q^{2}/2(1+q)}{\mathcal{F}}-\frac{q^{2}}{2(1+q)},$ (7)

where $\Phi_{\rm L1}$ is the potential at the L1 Lagrange point of the system.

We model the supernova explosion as a “thermal bomb” in a small region of the
grid,

$p(r)=\frac{3(\hat{\gamma}-1)E}{4\pi\delta r}\times\mathcal{H}(r-\delta r),$
(8)

where $\delta r$ is the radial extent of the overpressure, and
$\mathcal{H}(r)$ is the Heaviside step function. Since we are constrained
geometrically by the detailed evolution of the binary, we are only afforded
the liberty of varying $E/M$ and the Roche lobe filling factor, $\mathcal{F}$.
Results are presented with $2048$ radial zones per decade on a
logarithmically-spaced spherical-polar mesh with $\Delta
r=r\sin(\pi/2)\Delta\phi$. The radial extent is $r_{\rm min}=10^{-3}R_{\odot}$
to $r_{\rm max}=10R_{\odot}$ and the azimuthal extent covers half the sphere,
$\phi\in[0,\pi]$. This results in a grid of 8192 radial zones by 2795 angular
zones. In other words, the explosion is resolved in 361 radial zones, and the
star is resolved within 7644 radial zones. We set $\delta r=1.5r_{\rm min}$.

##  3 Results

Figure 1: The dynamical evolution of a $10^{51}\,\mathrm{erg}$ explosion
within a Roche-lobe filling stripped star ($\mathcal{F}=1$, shown as the grey
dashed line with volume-equivalent radius $\sim 3.8\times
10^{11}\,\mathrm{cm}$). Logarithmic pressure is shown on the right half and
four-velocity on the left half. In panel (a), the shock wave evolves quasi-
spherically while still inside the stellar mantle. In panel (b), the supernova
shock has broken out of the stellar poles first and begins wrapping around the
star towards collision points at the Roche-lobe nose and the back of the
progenitor. In panel (c), the lobe shocks collide and force material along the
center-of-mass axis of the binary, forming two relativistic bullets directed
towards and away from the companion with masses $6\times 10^{27}\,\mathrm{g}$
and $5\times 10^{25}\,\mathrm{g}$, of order the masses of the Earth and Moon,
respectively.

Figure 1 depicts the shock wave evolution from a $10^{51}\,\mathrm{erg}$
explosion set off inside the $3M_{\odot}$ progenitor. The first snapshot of
the explosion shows that at first, the supernova shock wave propagates outward
quasi-spherically while still interior to the stellar mantle in panel (a).
This is because the iso-density contours near the stellar core are quasi-
spherical due to gravity, so there is no preferential density gradient in
either direction that would initially accelerate and deform the supernova
shock (Sakurai, 1960; Matzner & McKee, 1999). However, because of the geometry
of the Roche lobe filling star, the shock waves will break out at the stellar
poles first, i.e., at $\theta=0$ and _not_ $\phi=0$, since the ratio of the
distance to one of the stellar poles to the distance to the back of the star
is $R_{z}/R_{\rm back}\approx 0.85$. The result of this asynchronous breakout
is shown in panel (b) of Figure 1 where the shock continues its outward
expansion, but the oblique breakout causes the shocks to wrap around towards
the nose of the Roche lobe as well as the back of the progenitor. The wrap-
around effect is caused by the naked blast wave realizing a lateral pressure
gradient in the $\hat{\phi}$ direction resulting in the blast wave forming a
pattern wave along the stellar surface, pointing towards impending collision
points along the center-of-mass axis of the binary (see e.g., Matzner et al.,
2013; Afsariardchi & Matzner, 2018; Scully et al., 2023, and references
therein). Finally, in panel (c) we see the aftermath of the lobe shocks
colliding toroidally at a point located along the center-of-mass line of the
binary system. The shocked lobes coalesce at this point, forming a hot,
pressurized blob of ejecta whose secondary acceleration is caused by the steep
pressure gradient directed along the binary center-of-mass axis. This
ultimately results in a violent spray of ejecta that form collimated
relativistic bullets pointing both away and towards the companion. The front
bullet has a mass of $6\times 10^{27}\,\mathrm{g}$ and the back bullet has a
mass of $5\times 10^{25}\,\mathrm{g}$, which are of order the masses of the
Earth and Moon, respectively.

To quantify some of the observable properties of the relativistic ejecta, a
useful measure is the cumulative kinetic energy distribution above some four-
velocity, $E_{k}(>\Gamma\beta)$. Figure 2 plots $E_{k}(>\Gamma\beta)$ for all
geometries and explosion energies considered in this work. There is non-
monotonic behavior in the Roche lobe stars where the geometries
$\mathcal{F}\geq 0.98$ extend the high-velocity tail beyond the ballistic
velocities shown for their spherical counterparts. The $\mathcal{F}=0.95$
stars, however, produce no material in such high-velocity ejecta and are
slower overall than the spherical counterparts. Although the $\mathcal{F}\geq
0.98$ explosions extend the high-velocity tail of the kinetic energy
distribution by factors of only a few, the angular distribution of this fast
tail shows a clear split in the azimuthal angle.

Figure 3 plots, for the $\mathcal{F}=1$ progenitors only, the cumulative
kinetic energy distribution per unit solid angle for the
$10^{51}\,\mathrm{erg}$ and $10^{52}\,\mathrm{erg}$ cases at various velocity
cutoffs. In panel (a), to demonstrate the global nature of the ejecta, we plot
the $\Gamma\beta=0$ cutoff and can recover the $10^{51}\,\mathrm{erg}$
explosion energy if multiplied by the full solid angle $4\pi$. However, when
we discriminate towards the mildly relativistic ($\Gamma\beta\simeq 0.5$) and
relativistic ($\Gamma\beta\simeq 1$) ejecta, we see energy is largely focused
along the binary axis for both explosion energies considered. For example, let
$dE_{k}/d\Omega\equiv E_{\Omega}$. In the $10^{51}\,\mathrm{erg}$ case we find
$E_{\Omega}({>}\\{0.5,0.1\\})|_{\phi=0^{\circ}}=\\{1.5\times
10^{45}\,\mathrm{erg},5\times 10^{44}\,\mathrm{erg}\\}$ and
$E_{\Omega}({>}0.5)|_{\phi=180^{\circ}}=4\times 10^{44}\,\mathrm{erg}$. In the
$10^{52}\,\mathrm{erg}$ case, we find that most of the energy is just focused
towards the companion with values
$E_{\Omega}({>}\\{0.5,1\\})|_{\phi=0^{\circ}}=\\{10^{48}\,\mathrm{erg},5\times
10^{47}\,\mathrm{erg}\\}$.

The solid angle of the bullets at the different velocity cutoffs is from

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle 4\pi\frac{E_{\rm
iso}}{E}$ (9) $\displaystyle=$ $\displaystyle
4\pi\left(\frac{\int_{\Omega}\left(dE/d\Omega\right)^{2}d\Omega}{\int_{\Omega}(dE/d\Omega)d\Omega}\right)\left(\int_{\Omega}(dE/d\Omega)d\Omega\right)^{-1},$

where $E_{\rm iso}$ is the isotropic-equivalent energy weighted by the energy
per solid angle. With Equation 9, we compute
$E_{k}({>}\\{0.5,1\\})|_{\phi=0^{\circ}}=\\{2\times
10^{44}\,\mathrm{erg},5\times 10^{41}\,\mathrm{erg}\\}$, and
$E_{k}({>}0.5)|_{\phi=180^{\circ}}=6\times 10^{43}\,\mathrm{erg}$ as the total
kinetic energies pointed towards and away from the companion for the
$10^{51}\,\mathrm{erg}$ explosion. For the $10^{52}\,\mathrm{erg}$ explosion,
the total energies are $E_{k}({>}\\{0.5,1\\})|_{\phi=0^{\circ}}=\\{4\times
10^{45}\,\mathrm{erg},2\times 10^{45}\,\mathrm{erg}\\}$.

Lastly, Figure 4 shows the separation between $10^{51}\,\mathrm{erg}$
explosions in the canonical spherical progenitor and a star that has filled
its Roche lobe. At freeze out of the acceleration in both cases, the Roche
lobe filling star produces trans-relativistic ejecta that gets focused along
the binary separation axis while there is no such fast ejecta in the perfectly
spherical supernova explosion. This distinctive asymmetry in the velocity and
density profiles might lend itself to various hints at the supernova
signatures observed and are further discussed in the next sections. We remind
the reader that these results are idealized to 2D axisymmetry, which may
enhance the effect of the shock converging at the nose and backside of the
Roche-lobe. Exploration of the 3D nature of this problem will be needed and we
plan to undertake this in a future study.

Figure 2: The ballistic cumulative kinetic energy distribution for the
$10^{51}\,\mathrm{erg}$ and $10^{52}\,\mathrm{erg}$ explosions at
$162\,\mathrm{s}$ and $58\,\mathrm{s}$, respectively, for all geometries
considered. The progenitors with Roche-lobe filling factors $\mathcal{F}\geq
0.98$ extend the high-velocity tail of the ejecta by factors of a few when
compared with their spherical counterparts. In general, the $\mathcal{F}=0.95$
produces the least amount of relativistic material. Figure 3: _Upper_ : The
cumulative kinetic energy per unit solid angle for the $10^{51}\,\mathrm{erg}$
explosion as a function of azimuthal angle, $\phi$, in the star that has
filled its Roche lobe, i.e., $\mathcal{F}=1$. The inset marks the region aimed
at the companion where the fastest ejecta is focused, carrying total energy
$E_{k}({>}1)=\Omega_{\rm bullet}dE_{k}/d\Omega=5\times 10^{41}\,\mathrm{erg}$
towards the companion in a bullet of angular width $\sim 1^{\circ}$. There is
also a counter bullet pointed away from the companion carrying away total
energy $\sim 10^{44}\,\mathrm{erg}$ in ejecta with $\Gamma\beta>0.5$. _Lower_
: The cumulative kinetic energy per unit solid angle for the
$10^{52}\,\mathrm{erg}$ explosion where most of the focusing happens towards
the companion with $E_{k}({>}1)=\Omega_{\rm bullet}dE_{k}/d\Omega=2\times
10^{45}\,\mathrm{erg}$. Figure 4: Evolution of the $10^{51}\,\mathrm{erg}$
supernova a few minutes post explosion in the stripped star (a) and in the
same stripped star that has filled its Roche lobe (b). Lab frame density is
plotted on the right side and four-velocity is plotted on the left. The grey-
dashed curves mark the original outer surface of the $3M_{\odot}$ progenitor
with volume-equivalent radius $3.8\times 10^{11}\,\mathrm{cm}$. The outflow
velocity for the spherically symmetric star peaks at around $\Gamma\beta=0.5$
while the fastest ejecta is focused along the binary axis for the Roche lobe
filling star and has velocity $\Gamma\beta\gtrsim 1$.

##  4 Astrophysical Implications

###  4.1 A relativistic bullet aimed at secondary and implications for shock
breakout

The late-stage morphology of explosion in the Roche-lobe distorted star
implies that the companion lies along the path of one of the relativistic
bullets. For the companion star to dodge the relativistic bullet, its
sectorial orbit would need to sweep an angle greater than the sum of the solid
angle of the companion with that of the bullet which leads to the companion
radial constraint:

$\frac{R_{\rm C}}{d}<\tan\left(\Delta\phi_{\rm orb}-\Delta\phi_{\rm
bullet}\right),$ (10)

where $R_{\rm C}$ is the radius of the companion, $d$ is the binary
separation, $\Delta\phi_{\rm orb}\equiv\tan^{-1}(V_{\rm orb}/\beta_{\rm
bullet}c)$ is the sectorial width of the companion orbit in the time it takes
the relativistic bullet to reach the companion’s original position, $V_{\rm
orb}\approx(GM/d)^{1/2}$ is the orbital velocity of the companion, $\beta_{\rm
bullet}$ is the velocity of the ballistic bullet, and $\Delta\phi_{\rm
bullet}$ is the half-opening angle of the bullet. With the bullet moving at
the speed of light and having $\Delta\phi_{\rm bullet}\simeq 0.01$, Inequality
10 requires $V_{\rm orb}>0.01c$ in order for the argument in the tangent
function to be non-negative. Since typical orbital velocities in binary
systems don not exceed a few hundred kilometers per second, this means that
the bullet will hit any companion regardless of its type whether it is a
neutron star (NS), a black hole (BH), or a main sequence (MS) star. We reserve
discussion of compact companions for later sections since their predicted
rates are lower than those of MS companions (see e.g., Zapartas et al., 2017;
Renzo et al., 2019).

While the energy budget of the directed bullet is not enough to unbind typical
envelopes for massive companion stars, the deposition of energy could
potentially affect the companion’s evolution through induced mixing effects,
seismic activity, or the collision could release a burst of radiation visible
from the shadow cone (Kasen, 2010; Ogata et al., 2021). Furthermore, the
companion star will be shock heated, expand, and radiate, which will produce
late-time radiation signatures at the supernova location (see e.g., Maund et
al., 2004; Folatelli et al., 2014; Maund et al., 2016; Ryder et al., 2018;
Maund, 2019; Sun et al., 2020, for up-to-date observations of stripped-
envelope SNe with companion detections), but since the energy imparted by the
bullet will be small, the companion’s inflation timescale will be shorter
lived than values expected for spherically symmetric SN hitting their
companions. We also note that if the explosion energy is large enough, the
ejecta-ejecta collision along the binary axis could produce a flash of photons
before the SN light is seen as discussed in Scully et al. (2023), and the
combined energies of the bullets in the $10^{52}\,\mathrm{erg}$ explosion
simulated in this work are $\sim 10^{46}\,\mathrm{erg}$, roughly the released
energy of XT 080901 (Soderberg et al., 2008). Moreover, the stripped-envelope
progenitors studied in this work are possible candidates for rapidly evolving
SNe now better known as a subclass of fast blue optical transients (FBOTs;
Drout et al., 2014; Kleiser & Kasen, 2014; Kleiser et al., 2018b, a), though
the luminosity peak will shift since the breakout flash is partially obscured
by the obliquity of the breakout (Matzner et al., 2013). All of these together
— the ejecta-ejecta flash, short shock cooling timescale of the companion, and
non-standard light curve morphology — might provide hints towards
classification of asymmetric SNe in binaries versus their spherically
symmetric counterparts.

###  4.2 Polarization signatures

In principle, polarization measurements can be made if there is asymmetry in
the supernova progenitor itself or if the explosion engine was asymmetric at
the outset (Höflich, 1991; Yamada & Sato, 1991; Höflich et al., 1996). In the
case for the Roche lobe bullets discussed in this work, the asymmetric outflow
will likely give rise to measurable non-zero net polarization, assuming that
it is dominated by Thompson scattering on free electrons. To get a sense of
the lifetime of the asymmetry in the flow, we estimate the deceleration of the
backwards bullet based on results from our simulations and assumptions about
the stellar environment. The front bullet will be quickly destroyed by
interactions with the companion, but the backwards bullet will decelerate
after a time $\left(\frac{M_{\rm bullet}}{4\pi A}\right)/\beta c$, which, when
scaled for a few reference parameters from the $10^{51}\,\mathrm{erg}$
explosion, gives

$t_{\rm dec}=12\left(\frac{M_{\rm bullet}}{5\times
10^{25}\,\mathrm{g}}\right)\left(\frac{A_{*}}{10^{-3}}\right)^{-1}\left(\frac{\beta}{0.5}\right)^{-1}\,\mathrm{days}$
(11)

and we have chosen $A_{*}=10^{-3}$ such that asymmetry of the outflow lasts
$\sim$ week which is roughly the earliest cadence at which we can currently
observe SNe post explosion (see Wang & Wheeler, 2008; Patat, 2017, and
references therin). Equation 11 is two-pronged. That is to say, the stellar
environment surrounding the Roche lobe explosion would have to be quite
rarefied or the mass of the bullet (and therefore the explosion energy) would
have to be larger in order to measure polarization signatures for these types
of explosions. At least one of these conditions can be satisfied based on: (i)
the large upper limits placed on the explosion energies of core-collapse
events; (ii) the stellar environment surrounding stripped-envelope stars being
uncertain. Furthermore, if the asymmetry lasts much longer than we estimate,
it is possible that bullets in binaries could produce core-collapse SN
remnants with “ears” (Grichener & Soker, 2017).

###  4.3 Asymmetric engines

Although this work is in the context of spherically symmetric explosions in
aspherical progenitors, there lies a plethora of observational evidence that
the resultant NS star formed from core-collapse can have a natal kick from an
asymmetric explosion mechanism (see e.g., Hobbs et al., 2005a; Scheck et al.,
2006; Wongwathanarat et al., 2013; Janka, 2017). An asymmetric explosion in an
asymmetric progenitor will likely greatly change the distribution of energy
over solid angle shown in Figure 3 since the resultant bullet velocities
strongly depend on the time and the angle at which they collide along the
binary axis. The advantage of the point explosion is the exterior lobe shocks
colliding simultaneously with the interior SN shock that breaks out at the
Roche lobe nose. This advantage of simultaneous collisions depends on the
explosion energy and $\mathcal{F}$ as shown in Figure 2. When the explosion
itself is asymmetric, the dynamics of the collision outside of the Roche lobe
filling star will change in a non-trivial way and should be studied in a
future work.

###  4.4 Fast Radio Bursts

In the scenario in which the binary companion is a neutron star (NS), the
relativistic bullet aimed toward the companion will encounter the strong
magnetic field of the NS magnetosphere. In this case, as the plasma bullet
hits the magnetic field and is decelerated, some of its kinetic energy can be
emitted in the form of radio waves. The mildly relativistic bullet has energy
$\sim 10^{44}\,\mathrm{erg}$ and initial radius of $\sim 3\times
10^{7}\,\mathrm{cm}$. Assuming the bullet is stopped almost instantaneously by
the strong NS magnetic field, then the time for the reverse shock to traverse
the bullet is roughly $10^{-3}\,\mathrm{s}$, of order the duration of a fast
radio burst (FRBs; Lorimer et al., 2007). We use the plasma frequency, $\nu$,
as an estimate for the radio wave frequency,

$\displaystyle\nu=\frac{\omega_{\rm plasma}}{2\pi}$ $\displaystyle=$
$\displaystyle\left(\frac{n_{e}\Omega_{\rm B}e^{2}}{\pi m_{e}\Omega_{\rm
bullet}}\right)^{1/2}$ (12) $\displaystyle=$ $\displaystyle
c\left(\frac{r_{e}n_{e}}{\pi}\right)^{1/2}\frac{(R_{\rm SO}/d)}{\theta_{0}},$

where $n_{e}=\rho/m_{p}$ is the electron number density in the bullet, $m_{p}$
is the proton mass, $e$ is the elementary charge, $m_{e}$ the electron mass,
$r_{e}$ is the classical electron radius, and we’ve accounted for the ratio of
solid angles of the bullet with that of the magnetosphere of the NS where
$R_{\rm SO}$ is the stand-off radius with value:

$\displaystyle R_{\rm SO}$ $\displaystyle\sim$ $\displaystyle R_{\rm
NS}\left(\frac{B_{d}^{2}}{8\pi\rho u^{2}c^{2}}\right)^{1/6}$ (13)
$\displaystyle\approx$ $\displaystyle 3R_{\rm
NS}B_{d,12}^{1/3}\rho_{-1}^{-1/6}u_{1}^{-1/3},$

where $B_{d}$ is the surface dipole field of the NS and we follow the
convention that quantities $Q_{x}=Q/10^{x}$ with cgs units implied. From this,
we find

$\nu\simeq 3\times
10^{12}B_{d,12}^{1/3}\rho_{-1}^{1/3}d_{12}^{-1}u_{1}^{-1/3}\theta_{0,-2}^{-1}\,\mathrm{Hz}.$
(14)

This implies that the bullet meets the energy, duration, and frequency
requirement for an FRB if it crashes into an NS magnetosphere, assuming an
efficiency $10^{-5}-10^{-4}$ for converting the kinetic energy into the
isotropic radio emission of $10^{39}\,\mathrm{erg}-10^{40}\,\mathrm{erg}$
quoted for many FRBs (see Cordes & Chatterjee, 2019, for a detailed review of
the current state of FRB science). Based on recently calculated rates of SNe
in the universe ($\sim 1$ s-1 Li et al., 2011) and FRB rates ($\sim 0.1$ s-1
Thornton et al., 2013) being in rough agreement with the rate of NS companions
existing in close binaries with massive stars (Zapartas et al., 2017; Renzo et
al., 2019), we predict that supernova explosions in close binaries with a
magnetized NS companion are a viable source for non-repeating FRBs via this
mechanism. In this scenario, the FRB would be followed by the rising light
curve of a stripped-envelope supernova (SN Type Ib/c) several days later.

###  4.5 Relevance for the formation of Double NS and GW sources

Explosions of Roche-Lobe distorted stars are expected in the main binary
formation scenarios for double compact objects. This is especially true in the
main formation pathway for double neutron stars (e.g. Dewi & Pols, 2003b;
Tauris et al., 2017), where both progenitors may be significantly distorted.
Our findings show potentially observable signatures that could be used to
search for systems evolving through this formation channel, see for example De
et al. (2018).

The first explosion is typically coming from the star that was once the most
massive star in the system. Depending on the separation and the amount of
swelling of the progenitor, this may already be an example of a Roche-
distorted star with an main sequence companion. Observable consequences may
result from ejecta hitting the companion, which is most likely still a main
sequence star in this phase (Sect. 4.1) and the aspherical outflows which may
give rise to polarization signatures (Sect. 4.2).

Most systems unbind during the first explosion, due the supernova kick
resulting from asymmetries in the explosion (e.g. Hobbs et al., 2005b).
Although we show in this work that additional asymmetries may result from the
deformation of the star, we find that these will not play any significant role
for the momentum imparted onto the system and the newly born neutron star.

Systems that remain bound can evolve further as a binary system. The second
star now evolves to fill its Roche-lobe. It first loses its hydrogen-rich
envelope in a highly non conservative event during which the orbit shrinks.
The result is a compact helium star orbiting a neutron star in a tight orbit.
The helium star is expected to swell once its start burning helium in a shell
around its core (often referred to as case BB mass transfer). It will fill the
Roche lobe anew now feeding gas to the neutron star, a stage that may be
observable as an X-ray binary. Eventually the helium star will complete the
advanced burning stages it will explode too, likely while still filling its
Roche-lobe. At this stage we expect a heavily Roche-lobe distorted star with a
neutron star companion, possibly with the observational consequence of a fast
radio burst (Sect. 4.4) accompanied by a supernova.

##  5 Summary and Conclusions

This work has demonstrated the following: (1) point-like explosions in pear-
shaped stripped-envelope stars cause shocks to converge to a point, producing
relativistic bullets pointed away and towards their companions even for the
typically non-relativistic $10^{51}\,\mathrm{erg}$ explosions; (2) such
explosions only impart $\sim 10^{44}-10^{45}$ erg onto their companions,
drastically reducing the shock cooling timescales for main sequence companions
assumed for spherical SN in binaries; (3) these explosions can obscure the
breakout flash and therefore the luminosity peak as seen for many rapidly
evolving SNe; (4) point-like explosions in pear-shaped stripped-envelope stars
might be candidates for Fast Radio Bursts associated with a core-collape SNe
if the relativistic bullet interacts with the magnetic field of a neutron star
companion. In short, as the vampire companion star depletes the primary of its
envelope, at the end of the primary’s life it responds by blasting the vampire
star with a relativistic bullet.

This work was supported in part by NASA ATP grant 80NSSC22K0822. MD and SdM
thank Rüdiger Pakmor, Mathieu Renzo, Shazrene Mohamed, Nando Patat, Rob
Farmer, Stephen Justham, Eva Laplace, and Norbert Langer for useful
discussions and comments during the MPA-Kavli summer school 2023, and MD
acknowledges useful discussions with Andrei Gruzinov and support from LSSTC
Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining
Grant #1829740, the Brinson Foundation, and the Moore Foundation.

## References

  * Afsariardchi & Matzner (2018) Afsariardchi, N., & Matzner, C. D. 2018, ApJ, 856, 146, doi: 10.3847/1538-4357/aab3d4
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Cordes & Chatterjee (2019) Cordes, J. M., & Chatterjee, S. 2019, ARA&A, 57, 417, doi: 10.1146/annurev-astro-091918-104501
  * De et al. (2018) De, K., Kasliwal, M. M., Ofek, E. O., et al. 2018, Science, 362, 201, doi: 10.1126/science.aas8693
  * Dewi & Pols (2003a) Dewi, J. D. M., & Pols, O. R. 2003a, MNRAS, 344, 629, doi: 10.1046/j.1365-8711.2003.06844.x
  * Dewi & Pols (2003b) —. 2003b, MNRAS, 344, 629, doi: 10.1046/j.1365-8711.2003.06844.x
  * Dewi et al. (2002) Dewi, J. D. M., Pols, O. R., Savonije, G. J., & van den Heuvel, E. P. J. 2002, MNRAS, 331, 1027, doi: 10.1046/j.1365-8711.2002.05257.x
  * Divine (1965) Divine, N. 1965, ApJ, 142, 824, doi: 10.1086/148356
  * Drout et al. (2023) Drout, M. R., Götberg, Y., Ludwig, B. A., et al. 2023, arXiv e-prints, arXiv:2307.00061, doi: 10.48550/arXiv.2307.00061
  * Drout et al. (2014) Drout, M. R., Chornock, R., Soderberg, A. M., et al. 2014, ApJ, 794, 23, doi: 10.1088/0004-637X/794/1/23
  * DuPont (2023) DuPont, M. 2023, SIMBI: 3D relativistic gas dynamics code, Astrophysics Source Code Library, record ascl:2308.003.  http://ascl.net/2308.003
  * DuPont et al. (2022) DuPont, M., MacFadyen, A., & Zrake, J. 2022, ApJ, 931, L16, doi: 10.3847/2041-8213/ac6ded
  * Efstathiou (2000) Efstathiou, G. 2000, MNRAS, 317, 697, doi: 10.1046/j.1365-8711.2000.03665.x
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960
  * Farmer et al. (2021) Farmer, R., Laplace, E., de Mink, S. E., & Justham, S. 2021, ApJ, 923, 214, doi: 10.3847/1538-4357/ac2f44
  * Filippenko (1997) Filippenko, A. V. 1997, ARA&A, 35, 309, doi: 10.1146/annurev.astro.35.1.309
  * Folatelli et al. (2014) Folatelli, G., Bersten, M. C., Benvenuto, O. G., et al. 2014, ApJ, 793, L22, doi: 10.1088/2041-8205/793/2/L22
  * Fowler & Hoyle (1964) Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201, doi: 10.1086/190103
  * Götberg et al. (2017) Götberg, Y., de Mink, S. E., & Groh, J. H. 2017, A&A, 608, A11, doi: 10.1051/0004-6361/201730472
  * Götberg et al. (2023) Götberg, Y., Drout, M. R., Ji, A. P., et al. 2023, arXiv e-prints, arXiv:2307.00074, doi: 10.48550/arXiv.2307.00074
  * Grichener & Soker (2017) Grichener, A., & Soker, N. 2017, MNRAS, 468, 1226, doi: 10.1093/mnras/stx534
  * Habets (1986) Habets, G. M. H. J. 1986, A&A, 167, 61
  * Hobbs et al. (2005a) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005a, MNRAS, 360, 974, doi: 10.1111/j.1365-2966.2005.09087.x
  * Hobbs et al. (2005b) —. 2005b, MNRAS, 360, 974, doi: 10.1111/j.1365-2966.2005.09087.x
  * Höflich (1991) Höflich, P. 1991, A&A, 246, 481
  * Höflich et al. (1996) Höflich, P., Wheeler, J. C., Hines, D. C., & Trammell, S. R. 1996, ApJ, 459, 307, doi: 10.1086/176894
  * Janka (2017) Janka, H.-T. 2017, ApJ, 837, 84, doi: 10.3847/1538-4357/aa618e
  * Kasen (2010) Kasen, D. 2010, ApJ, 708, 1025, doi: 10.1088/0004-637X/708/2/1025
  * Kippenhahn & Weigert (1967) Kippenhahn, R., & Weigert, A. 1967, ZAp, 65, 251
  * Kleiser et al. (2018a) Kleiser, I., Fuller, J., & Kasen, D. 2018a, MNRAS, 481, L141, doi: 10.1093/mnrasl/sly180
  * Kleiser & Kasen (2014) Kleiser, I. K. W., & Kasen, D. 2014, MNRAS, 438, 318, doi: 10.1093/mnras/stt2191
  * Kleiser et al. (2018b) Kleiser, I. K. W., Kasen, D., & Duffell, P. C. 2018b, MNRAS, 475, 3152, doi: 10.1093/mnras/stx3321
  * Klencki et al. (2022) Klencki, J., Istrate, A., Nelemans, G., & Pols, O. 2022, A&A, 662, A56, doi: 10.1051/0004-6361/202142701
  * Kopal (1959) Kopal, Z. 1959, Close binary systems
  * Laplace et al. (2020) Laplace, E., Götberg, Y., de Mink, S. E., Justham, S., & Farmer, R. 2020, A&A, 637, A6, doi: 10.1051/0004-6361/201937300
  * Leahy & Leahy (2015) Leahy, D. A., & Leahy, J. C. 2015, Computational Astrophysics and Cosmology, 2, 4, doi: 10.1186/s40668-015-0008-8
  * Li et al. (2011) Li, W., Chornock, R., Leaman, J., et al. 2011, MNRAS, 412, 1473, doi: 10.1111/j.1365-2966.2011.18162.x
  * Lorimer et al. (2007) Lorimer, D. R., Bailes, M., McLaughlin, M. A., Narkevic, D. J., & Crawford, F. 2007, Science, 318, 777, doi: 10.1126/science.1147532
  * Matzner et al. (2013) Matzner, C. D., Levin, Y., & Ro, S. 2013, ApJ, 779, 60, doi: 10.1088/0004-637X/779/1/60
  * Matzner & McKee (1999) Matzner, C. D., & McKee, C. F. 1999, ApJ, 510, 379, doi: 10.1086/306571
  * Maund (2019) Maund, J. R. 2019, ApJ, 883, 86, doi: 10.3847/1538-4357/ab2386
  * Maund et al. (2016) Maund, J. R., Pastorello, A., Mattila, S., Itagaki, K., & Boles, T. 2016, ApJ, 833, 128, doi: 10.3847/1538-4357/833/2/128
  * Maund et al. (2004) Maund, J. R., Smartt, S. J., Kudritzki, R. P., Podsiadlowski, P., & Gilmore, G. F. 2004, Nature, 427, 129, doi: 10.1038/nature02161
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Ogata et al. (2021) Ogata, M., Hirai, R., & Hijikawa, K. 2021, MNRAS, 505, 2485, doi: 10.1093/mnras/stab1439
  * Patat (2017) Patat, F. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 1017, doi: 10.1007/978-3-319-21846-5_110
  * 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
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66, doi: 10.1051/0004-6361/201833297
  * Ryder et al. (2018) Ryder, S. D., Van Dyk, S. D., Fox, O. D., et al. 2018, ApJ, 856, 83, doi: 10.3847/1538-4357/aaaf1e
  * Sakurai (1960) Sakurai, A. 1960, Communications on Pure and Applied Mathematics, 13, 353, doi: https://doi.org/10.1002/cpa.3160130303
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Scheck et al. (2006) Scheck, L., Kifonidis, K., Janka, H. T., & Müller, E. 2006, A&A, 457, 963, doi: 10.1051/0004-6361:20064855
  * Scully et al. (2023) Scully, B., Matzner, C. D., & Yalinewich, A. 2023, MNRAS, 525, 1562, doi: 10.1093/mnras/stad2360
  * Soderberg et al. (2008) Soderberg, A. M., Berger, E., Page, K. L., et al. 2008, Nature, 453, 469, doi: 10.1038/nature06997
  * Sravan et al. (2019) Sravan, N., Marchant, P., & Kalogera, V. 2019, ApJ, 885, 130, doi: 10.3847/1538-4357/ab4ad7
  * Sun et al. (2020) Sun, N.-C., Maund, J. R., Hirai, R., Crowther, P. A., & Podsiadlowski, P. 2020, MNRAS, 491, 6000, doi: 10.1093/mnras/stz3431
  * Tauris et al. (2015) Tauris, T. M., Langer, N., & Podsiadlowski, P. 2015, MNRAS, 451, 2123, doi: 10.1093/mnras/stv990
  * 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
  * Thornton et al. (2013) Thornton, D., Stappers, B., Bailes, M., et al. 2013, Science, 341, 53, doi: 10.1126/science.1236789
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004, doi: 10.21105/joss.02004
  * Wang & Wheeler (2008) Wang, L., & Wheeler, J. C. 2008, ARA&A, 46, 433, doi: 10.1146/annurev.astro.46.060407.145139
  * Wongwathanarat et al. (2013) Wongwathanarat, A., Janka, H. T., & Müller, E. 2013, A&A, 552, A126, doi: 10.1051/0004-6361/201220636
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015, doi: 10.1103/RevModPhys.74.1015
  * Yamada & Sato (1991) Yamada, S., & Sato, K. 1991, ApJ, 382, 594, doi: 10.1086/170746
  * Yoon et al. (2010) Yoon, S. C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940, doi: 10.1088/0004-637X/725/1/940
  * Zapartas et al. (2017) Zapartas, E., de Mink, S. E., Van Dyk, S. D., et al. 2017, ApJ, 842, 125, doi: 10.3847/1538-4357/aa7467


# Gel-like granular materials with high durability and high deformability

Honoka Fujio, Hikari Yokota, Marie Tani and Rei Kurita  Department of Physics,
Tokyo Metropolitan University, 1-1 Minamioosawa, Hachiouji-shi, Tokyo
192-0397, Japan

###### Abstract

Building materials such as concretes and mortar are formed by solidifying
granular slabs. Such materials are often fractured by giant forces such as
during earthquakes, leading to the collapse of structures and potentially
casualties. One avenue of enquiry to prevent cracking would be to realize a
material that can maintain a stable shape without being solidified. Here, we
focus on sand grains coated with silicone oil, experimentally investigating
the Young’s modulus of a granular slab where ordinary grains and the coated
grains are mixed in a mixing ratio $\alpha$. It is found that the Young’s
modulus increases rapidly at $\alpha\geq 0.6$. We use numerical simulation to
show that this sudden increase in the Young’s modulus is caused by a rigidity
percolation transition. Furthermore, we are able to show that granular slabs
containing coated sand have outstanding deformability without collapsing under
large external stress. We believe this may lead to the development of granular
materials that are rigid under usual pressures but deformable under more
extreme conditions, such as during seismic activity.

Granular systems are ubiquitously observed in daily life. Examples include
sandboxes, dunes, flour, beads, and pharmaceutical powders. Despite individual
particles obeying classical equations of motion, granular systems are known to
exhibit complex macroscopic behavior. They can behave like a fluid, as in an
hourglass, or like a solid, when they form arches or domes, clogging hoppers
and pipes Duran (2012); Brown and Richards (2016); Nagel (1992); Jaeger and
Nagel (1992); Guyon and Troadec (1994). Understanding the behavior of granular
systems is an important challenge for agriculture, construction,
pharmaceuticals, and preparing for natural disasters such as landslides and
avalanches.

It is known that forces inside dry granular systems do not propagate uniformly
but in a chain-like manner. This force chain is intricately linked to their
macroscopic behavior Drescher and De Jong (1972). Although dry granular
materials are well studied, it should be noted that granular materials found
in nature and our daily lives often contain liquid. With the addition of even
a small amount of liquid, the mechanical properties of granular materials are
dramatically altered. The liquid bridges grains (liquid adhesion or capillary
bridging), and the capillary force causes an effective attractive interaction
between the particles Mitarai and Nori (2006); Samadani and Kudrolli (2000);
Scheel  _et al._ (2008). Therefore, it has greater rigidity than dry granular
materials. Concretes, for example, are made by solidifying the liquid bridged
state, and are found in many everyday applications, such as constructions and
roads. Such materials do not have a relaxation time scale since they are
solidified, and consequently cannot relax their internal structure. When
gigantic perturbing forces such as earthquakes are applied, they can undergo
brittle fracture, resulting in human casualties. Therefore, there is a need
for novel granular materials that not only have high durability but also large
deformability.

Sand coated with silicone oil (coated sand) is a product designed as a
children’s toy Modell and Thuresson (2011). The silicone oil is bonded
directly to each grain and does not peel off. The coating allows the oil to
persist stably around the sand. When the oil layers come into contact with
each other, liquid adhesion occurs, causing an attractive interaction between
coated particles Tani  _et al._ (2021). Here we note that this attractive
interaction remains for a long time, while the capillary bridging cannot
retain due to evaporation or drainage. Experiments where this coated sand was
mixed with ordinary grains have been reported as a model system for wet
granular materials Tani  _et al._ (2021). In such systems, microscopy
observations have confirmed that no stable liquid adhesion occurs between
coated and uncoated sand, no matter how many times or how long they are mixed.
It was reported that the packing fraction of random loose packing, the angle
of repose, and the amount of residual material when the mixture is sieved
change significantly when the mixing ratio of the coated sand to the total
mass $\alpha$ is over 0.23 Tani  _et al._ (2021). This is close to the
connectivity percolation value of $\alpha\sim 0.25$ Stauffer and Aharony
(1994), suggesting that the change in physical properties may be due to
percolation of clusters formed by coated sand.

Previous experiments focused on physical properties with a mixing ratio of
$\alpha<0.3$, where the granular slab cannot retain its shape without a
container Tani  _et al._ (2021). It does not become a block and cannot be used
as a granular material. In this study, we focus on the mechanical properties
of randomly close packed (RCP) blocks of material as a possibility for actual
material development. Here we investigate the $\alpha$ dependence of the
Young’s modulus and deformability of a granular block. We also seek to clarify
the origins of large Young’s modulus and large deformability using numerical
simulations.

For ordinary grains, we used Tohoku Silica Sand No. 5 sand (Tohoku Silica Sand
Co., Ltd.). These grains have a density of 2.62 g/cm3, a short diameter of
0.08-0.35 mm, and a long diameter of 0.23-0.55 mm. Kinetic Sand (Rangs Japan)
was used for silicone oil coated sand. The density is 2.26 g/cm3, the short
diameter is 0.09-0.29 mm, and the long diameter is 0.16-0.41 mm. For each mass
ratio $\alpha$, we measured a mass of each, spread them on the same metal vat
and mixed them well by hand until uniform. Microscopic observation of several
locations confirmed that the mixture was nearly uniform.

We prepared the random close packed (RCP) slabs in the following way. We
gently put the sand into a cylindrical container with a diameter of 74.0 mm
and height of 42.0 mm, or a rectangular container with a 67.8 mm $\times$ 67.8
mm base and a height of 42.5 mm. The rate at which the sand was poured was
about 10 g/s. RCP states are generated by compacting mixtures down by external
pressure ($\sim$ 5 $\times$ 104 Pa) using a plate whose area is almost the
same as the inside of the container. We then poured more sand over the top
until the container is full. We repeated this process around 25 times. The
packing fraction of the granular slabs become roughly uniform. This state can
be regarded as being RCP.

We flipped the container over and extracted the RCP granular slab from the
container (see Fig. 1). Next, external pressure was applied by using a slider
(LT-S300, Thorlabs) equipped with a large plate (Fig. 2(a)). The plate is
designed so that it always stays horizontal. The slider can move with a
constant velocity $v$ with a minimum $v$ of 0.001 mm/s. The force was measured
using a force gauge (DST-50N, IMADA, Japan) during the slider movement. The
resolution of the force is 0.001 N. We measured this 10 times using the same
experimental conditions; the error bars correspond to the maximum and minimum
found during the experimental trials.

We performed overdamped numerical simulations to investigate the internal
structures of granular slabs. To simplify, we assume grains to be spherical
with a diameter $d_{i}$ for particle $i$. We set the size distribution to be
Gaussian with a polydispersity of 20 %. A two-body potential $U(r_{ij})$
between particles $i$ and $j$ is applied such that
$U(r_{ij})=\frac{k}{2}[r_{ij}-(d_{i}+d_{j})/2]^{2}$ if
$(d_{i}+d_{j})/2>r_{ij}$, where $r_{ij}$ is the distance between particles $i$
and $j$ and $k$ = 1. Otherwise, $U(r_{ij})=0$. Firstly, we prepared a RCP
state of non-coated grains using a method reported in Ref. Xu  _et al._
(2005); Clarke and Wiley (1987). The particles are overpacked until the
packing fraction is 0.645, and then shrunk until the total potential is below
10-3. Next, we replace a certain number of non-coated sand grains with coated
grains until the number ratio of coated grains becomes $\alpha$. Coated grains
$i$ and $j$ interact via an attractive potential
$U(r_{ij})=\frac{k_{att}}{2}[r_{ij}-(d_{i}+d_{j})/2]^{2}$ when
$(d_{i}+d_{j})/2<r_{ij}<r_{c}$, where $r_{c}$ is a threshold for the
attractive interaction. We set $k_{att}$ = 0.1 and
$r_{c}=1.05(d_{i}+d_{j})/2$; we also confirm that our results are not
sensitive to the threshold value. We measured 3 different particle
configurations and 5 different ways to replace uncoated grains with coated
grains for each data point.

Firstly, we examined the stability of cylindrical RCP samples under their own
weight. Figure 1 shows the appearance of the samples with $\alpha$ = 0.0, 0.3,
and 0.5, 60 minutes after it was removed from a cylindrical container. When
$\alpha\leq$ 0.25, the sample collapses the moment it is removed from the
container and forms a mound [Fig. 1(a)]. The angle of repose is 20∘ $\sim$
25∘, which is almost the same as that of an ordinary sand pile Tani  _et al._
(2021); Brown and Richards (2016). On the other hand, a slab with
$0.3\leq\alpha\leq 0.4$ only partially collapses as shown in Fig. 1(b), almost
maintaining its shape. Finally, when the ratio of coated sand increases to
$\alpha\geq 0.5$, the granular slab maintains its shape completely, as shown
in Fig. 1(c). The slab keeps its shape perfectly and does not collapse after a
long time.

Figure 1: Appearance of samples 60 minutes after they were removed from a
cylindrical container, where (a) $\alpha$ = 0.0, (b) 0.3, and (c) 0.5. The
scale bar corresponds to 50 mm.

Next, we investigate response to large external forces. As shown in the
experimental schematic in Fig. 2(a), samples were compressed from the top
surface at a constant speed by a plate, and the normal stress $\sigma$ was
measured. When $\alpha\leq 0.3$, the sample collapsed immediately and $\sigma$
could not be measured. When $\alpha\geq 0.4$, they did not collapse under
small deformations. Figure 2(b) shows the dependence of the normal stress
$\sigma$ on strain $\varepsilon$ for $\alpha$ = 0.4, 0.6, 0.8 and 1.0 at a
compression speed of $v$ = 0.01 mm/s. When the strain is small, the normal
stress increases linearly; as the strain is increased further, the stress
increases rapidly. However, we see that the $\sigma$-$\varepsilon$ curve is
qualitatively different for $\alpha$ = 0.4 and $\alpha\geq 0.6$. Here, we
measured the slope of the linear region under small strains and estimated the
Young’s modulus of the grain slab, $E=\sigma/\varepsilon$, from the slope.
Figure 2(c) shows the relationship between $\alpha$ and $E$. The Young’s
modulus changes rapidly for $0.5<\alpha<0.6$ and then remains almost constant
for $\alpha\geq 0.6$.

The same experiment was conducted at different compression speeds. As shown in
Fig. 2(c), the Young’s modulus increases sharply at $\alpha\sim 0.6$ for all
$v$. The means of the Young’s modulus for 0.4 $\leq\alpha\leq$ 0.5 and at 0.6
$\leq\alpha\leq$ 1 are calculated as shown in Fig. 2(c) and labeled $E_{1}$
and $E_{2}$, respectively for each $v$. $E_{1}$ and $E_{2}$ increased sharply
by several times at a transition compression rate $v_{t}\sim 0.2$ mm/s as
shown in Fig. 2(d). This suggests the existence of a relaxation time scale for
the internal structure. The characteristic time may be estimated by
$\tau\sim\langle d\rangle/v_{t}~{}\sim$ 1 s, where $\langle d\rangle$ is the
mean size of the grains. A possible explanation for this time scale can be the
lifetime of the attractive interaction. The time scale over which this liquid
bridge disconnects is about 1 s for our coated grains. Another possibility is
structural relaxation, where particles in contact rearrange in the same way as
$T1$ events in foams Yanagisawa and Kurita (2021). We are also working on a
detailed study of the origin of this time scale and will report these results
in the near future.

Figure 2: (a) Schematic for the experimental setup, viewed from the side. (b)
$\sigma$-$\varepsilon$ curve. When the strain is small, the normal stress
increases linearly (dashed line) and the Young’s modulus of the grain slab $E$
can be estimated as $E=\sigma/\varepsilon$ from the slope. (c) Young’s modulus
$E$ as a function of $\alpha$. The Young’s modulus increases sharply at
$\alpha\sim 0.6$ for all compression rates $v$. The mean values of the Young’s
modulus over 0.4 $\leq\alpha\leq$ 0.5 and 0.6 $\leq\alpha\leq$ 1 are labeled
$E_{1}$ and $E_{2}$, respectively for each $v$. The solid lines represent
$E_{1}$ and $E_{2}$ for $v$ = 0.5 mm/s. (d) The Young’s modulus as a function
of the compression rate $v$. $E_{1}$ and $E_{2}$ increase sharply by
approximately 5 times at a transition compression rate of $v_{t}=0.2$ mm/s.

Here, we note the importance of the granular slab having a structural
relaxation in the context of construction materials. Normal concretes undergo
brittle fracture when they are deeply compressed or a crack emerges inside.
However, granular slabs with $\alpha\geq 0.6$ can be deformed without cracking
even when $\varepsilon=0.15$, as shown in Fig. 3(a1)-(a3); it is found that
the width at the center of the granular slab expanded more than at the base.
This means that the structure may have changed inside the slab, something
which cannot occur when the slab is solidified. We also check the
deformability for the shear deformation. We sandwiched the granular slab by
two acrylic plates and then we sheared the granular slab by tilting the
acrylic plates. The granular slabs with $\alpha\geq 0.6$ can be deformed even
when $\varepsilon=0.15$, as shown in Fig. 3(b1)-(b3). The granular slabs with
$\varepsilon=0.15$ is stable at least for an hour.

Figure 3: Snapshots of a granular slab with $\alpha=0.6$ under different
strains, $\varepsilon=$ (a1) 0.0, (a2) 0.10, and (a3) 0.15. Dotted and dashed
lines show the positions of the top surface at $\varepsilon$ = 0.0, and 0.1,
respectively. Snapshots of a rectangular granular slab with $\alpha=0.6$ under
different shear strains $\varepsilon=$ (b1) 0.0, (b2) 0.10, and (b3) 0.15.
Dotted vertical lines shows the perpendicular to the substrate and dash lines
show the tilted edge of the granular slab. The granular slab can be deformed
without cracking even though $\varepsilon=0.15$ and the granular slab is
stable at least for an hour. The scale bar corresponds to 50 mm.

Next, to clarify the origin of the change of the Young’s modulus, we compute
the inner structures using a numerical simulation of packings spheres which
model our mixture; details given in the Methods section. Firstly, we focus on
connections between coated sand grains. Figure 4 shows the spatial
distribution of clusters composed of coated sand grains at (a) $\alpha$ = 0.2
and (b) 0.3. The colors represent the size of each cluster. It is found that
many small clusters are scattered when $\alpha$ = 0.2, while there is one
large cluster when $\alpha$ = 0.3. We measured the particle number in each
connected cluster $N_{c}$ and defined $N_{c}^{max}$ as the number of particles
in the largest cluster. Figure 4(c) shows the ratio of $N_{c}^{max}$ to the
total number of coated sand particles $\beta_{c}$, that is,
$\beta_{c}=N_{c}^{max}/N\alpha$. It is found that $\beta_{c}$ increases
significantly at $\alpha=0.3$. As $\alpha$ increases further, $\beta_{c}$
approaches 1; this means that almost all coated sand grains are included in
one large cluster. We also confirm that the large cluster percolates the
system. Furthermore, the contact number $Z$ in the maximum cluster is $Z$ =
2.77 at $\alpha=0.3$. It has been reported previously that percolated
structures are rigid when $Z\geq 2.4$ Kuo and Gupta (1995); Gallegos  _et al._
(2021). Thus, the connectivity percolation transition corresponds to the
shape-retaining property of the granular slab shown in Fig. 1. This is also
consistent with the correspondence of the connectivity percolation with
changes in the packing fraction of random loose packing, behavior under
sieving, and the angle of repose as reported in Ref. Tani  _et al._ (2021).

Figure 4: Spatial distribution of clusters composed of connected coated sand
grains when (a) $\alpha$ = 0.2 and (b) 0.3. The colors correspond to the size
of each cluster. It is found that many small clusters are scattered over space
when $\alpha$ = 0.2, while one large cluster percolates the system at $\alpha$
= 0.3. (c) $\beta_{c}$ as a function of $\alpha$. The maximum cluster size
increases significantly at $\alpha=0.3$. The error bars correspond to the
standard deviations.

The connectivity percolation should increase the stiffness of the granular
slab, but it is also known that connectivity is not enough to realize
mechanical rigidity Maxwell (1864); Guyon  _et al._ (1990); Aharonov and
Sparks (1999); Lois  _et al._ (2008); Yanagishima  _et al._ (2021). For
example, when a particle makes contact with two particles, a linear cluster
can be formed. Such a cluster is weak against large stresses and can be easily
deformed. Meanwhile, Maxwell’s condition stipulates that connections become
rigid when the total contact number $z$ reaches $6n$-12, where $n$ is the
particle number Maxwell (1864). Therefore, the contact number for each
particles is almost 6 since $n$ is enough large. We defined a rigid connection
if the contact number is over 6. We apply this to the structures generated in
our simulation and investigate the spatial distribution of rigid connections.
Figure 5 shows the distribution of rigid clusters at (a) $\alpha$ = 0.5 and
(b) 0.6. Now, colors correspond to the size of rigid clusters. It is found
that small rigid clusters are scattered over space at $\alpha$ = 0.5, even
though the clusters as judged by connections are percolated. This means that
the connective clusters are deformable when a large external force is applied.
Meanwhile, when $\alpha$ = 0.6, it is clear that there is one large rigid
cluster. We measured the number of particles in each rigid cluster $N_{r}$ and
define $N_{r}^{max}$ as the particle number in the largest cluster. Figure
5(c) shows the ratio of $N_{r}^{max}$ to the total number of coated particles
$\beta_{r}$, that is, $\beta_{r}=N_{r}^{max}/N\alpha$. It is found that
$\beta_{r}$ becomes non-zero at $\alpha=0.6$; the value $\beta_{r}$ = 0.3 is
larger than that required for percolation; we also directly confirm that the
rigid cluster percolates the system by checking the coordinates of the
particles. This rigidity percolation is the reason why the Young’s modulus
suddenly increases for $\alpha\geq 0.6$.

Figure 5: Spatial distribution of rigid clusters where the contact number is
over 6 when (a) $\alpha$ = 0.5, (b) $\alpha$ = 0.6. The color represents the
size of each rigid cluster. It is found that many small rigid clusters are
scatter over space at $\alpha$ = 0.5, while one large rigid cluster percolates
the system at $\alpha$ = 0.6. (c) $\beta_{r}$ as a function of $\alpha$. Note
that $\beta_{r}\sim$ 0.3 at $\alpha=0.6$; this is larger than the value
required for percolation. The error bars correspond to the standard
deviations.

Here, the mechanical properties of this mixed granular slab strongly depend on
the connection between coated sand particles via an attractive interaction.
This has a close analogy with gels, where polymer networks may percolate
Okumura and Ito (2001); Sakai  _et al._ (2008). Meanwhile, we note that
granular systems are generally related to jammed system. Thus, this mixed
granular system may have the physical properties of both gels and granular
systems, such as large rigidity and large deformability without brittle
fracture. This may mean that new functional granular materials may be produced
by simply changing the coating of a certain proportion of grains.

We investigated a unique granular slab composed of a mixture of sand grains
coated and not coated with silicone oil. We measured the mechanical properties
of RCP granular slabs as a function of the mass ratio of coated sand grains
$\alpha$. It is found that granular slabs can retain their own shape when
$\alpha\geq 0.3$. We also found that the Young’s modulus suddenly increases at
$\alpha$ = 0.6. From the results of numerical simulations, the coated grains
not only percolate via connectivity at $\alpha\geq 0.3$, but also undergo
rigidity percolation when $\alpha\geq 0.6$. The former percolation is related
with shape retention, while the latter is correlated with the increased
Young’s modulus. Furthermore, the Young’s modulus depends on a compression
velocity, suggesting the granular slab has a characteristic relaxation time
scale. Thus, we can expect viscoelastic behavior in our granular slabs, paving
the way for new functional granular materials simply by changing the coating
of a proportion of grains.

M. T. was supported by JSPS KAKENHI (20K14431) and R. K. was supported by JSPS
KAKENHI (17H02945 and 20H01874). Correspondence and requests for materials
should be addressed to R. K. (kurita@tmu.ac.jp).

## References

  * Duran (2012) J. Duran,  _Sands, powders, and grains: an introduction to the physics of granular materials_ (Springer Science & Business Media, 2012).
  * Brown and Richards (2016) R. L. Brown and J. C. Richards,  _Principles of powder mechanics: essays on the packing and flow of powders and bulk solids_ , Vol. 10 (Elsevier, 2016).
  * Nagel (1992) S. R. Nagel, Reviews of Modern Physics 64, 321 (1992).
  * Jaeger and Nagel (1992) H. M. Jaeger and S. R. Nagel, Science 255, 1523 (1992).
  * Guyon and Troadec (1994) É. Guyon and J.-P. Troadec,  _Du sac de billes au tas de sable_ (Odile Jacob, 1994).
  * Drescher and De Jong (1972) A. Drescher and G. D. J. De Jong, Journal of the Mechanics and Physics of Solids 20, 337 (1972).
  * Mitarai and Nori (2006) N. Mitarai and F. Nori, Advances in Physics 55, 1 (2006).
  * Samadani and Kudrolli (2000) A. Samadani and A. Kudrolli, Physical Review Letters 85, 5102 (2000).
  * Scheel  _et al._ (2008) M. Scheel, R. Seemann, M. Brinkmann, M. Di Michiel, A. Sheppard, B. Breidenbach, and S. Herminghaus, Nature Materials 7, 189 (2008).
  * Modell and Thuresson (2011) J. Modell and S. Thuresson, “Material composition and method for its manufacture,” (2011), uS Patent 7,862,894.
  * Tani  _et al._ (2021) M. Tani, H. Fujio, and R. Kurita, J. Phys. Soc. Jpn. 90, 033801 (2021).
  * Stauffer and Aharony (1994) D. Stauffer and A. Aharony,  _Introduction To Percolation Theory_ (Taylor & Francis, 1994).
  * Xu  _et al._ (2005) N. Xu, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev. E 71, 061306 (2005).
  * Clarke and Wiley (1987) A. S. Clarke and J. D. Wiley, Phys. Rev. B 35, 7350 (1987).
  * Yanagisawa and Kurita (2021) N. Yanagisawa and R. Kurita, Sci. Rep. 11, 2786 (2021).
  * Kuo and Gupta (1995) C.-H. Kuo and P. K. Gupta, Acta Metallurgica et Materialia 43, 397 (1995).
  * Gallegos  _et al._ (2021) J. A. S. Gallegos, R. Perdomo-Pérez, N. E. Valadez-Pérez, and R. Castañeda Priego, Phys. Rev. E 104, 064606 (2021).
  * Maxwell (1864) J. C. Maxwell, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 27, 294 (1864).
  * Guyon  _et al._ (1990) E. Guyon, S. Roux, A. Hansen, D. Bideau, J.-P. Troadec, and H. Crapo, Reports on Progress in Physics 53, 373 (1990).
  * Aharonov and Sparks (1999) E. Aharonov and D. Sparks, Physical Review E 60, 6890 (1999).
  * Lois  _et al._ (2008) G. Lois, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev. Lett. 100, 028001 (2008).
  * Yanagishima  _et al._ (2021) T. Yanagishima, J. Russo, R. P. Dullens, and H. Tanaka, Phys. Rev. Lett. 127, 215501 (2021).
  * Okumura and Ito (2001) Y. Okumura and K. Ito, Advanced materials 13, 485 (2001).
  * Sakai  _et al._ (2008) T. Sakai, T. Matsunaga, Y. Yamamoto, C. Ito, R. Yoshida, S. Suzuki, N. Sasaki, M. Shibayama, and U.-i. Chung, Macromolecules 41, 5379 (2008).