We present unstratified 3D MHD simulations of an accretion disk with a boundary layer (BL) that have a duration $\sim 1000$ orbital periods at the inner radius of the accretion disk. We find the surprising result that angular momentum piles up in the boundary layer, which results in a rapidly rotating belt of accreted material at the surface of the star. The angular momentum stored in this belt increases monotonically in time, which implies that angular momentum transport mechanisms in the BL are inefficient and do not couple the accretion disk to the star. This is in spite of the fact that magnetic fields are advected into the BL from the disk and supersonic shear instabilities in the BL excite acoustic waves. In our simulations, these waves only carry a small fraction ($\sim 10 \%$) of the angular momentum required for steady state accretion. Using analytical theory and 2D viscous simulations in the $R-\phi$ plane, we derive an analytical criterion for belt formation to occur in the BL in terms of the ratio of the viscosity in the accretion disk to the viscosity in the BL. Our MHD simulations have a dimensionless viscosity ($\alpha$) in the BL that is at least a factor of $\sim 100$ smaller than that in the disk. We discuss the implications of these results for BL dynamics and emission.

PICsar2D is a 2.5D relativistic, electromagnetic, particle in cell code designed for studying the pulsar magnetosphere. The source code and a suite of Python analysis routines can be downloaded from "https://github.com/astromb/PICsar2D.git". Additionally, the repository includes a step-by-step tutorial for compiling the code, running it, and analyzing the output. This article is devoted to several new algorithmic advances and numerical experiments. These include a new pair injection prescription at the pulsar surface, a comparison of different pair injection techniques, a discussion of particle trajectories near the pulsar Y-line, and performance optimization of the code.

We present an instability for exciting incompressible modes (e.g. gravity or Rossby modes) at the surface of a star accreting through a boundary layer. The instability excites a stellar mode by sourcing an acoustic wave in the disk at the boundary layer, which carries a flux of energy and angular momentum with the opposite sign as the energy and angular momentum density of the stellar mode. We call this instability the acoustic CFS instability, because of the direct analogy to the Chandrasekhar-Friedman-Schutz instability for exciting modes on a rotating star by emission of energy in the form of gravitational waves. However, the acoustic CFS instability differs from its gravitational wave counterpart in that the fluid medium in which the acoustic wave propagates (i.e.\ the accretion disk) typically rotates faster than the star in which the incompressible mode is sourced. For this reason, the instability can operate even for a non-rotating star in the presence of an accretion disk. We discuss applications of our results to high-frequency quasi-periodic oscillations in accreting black hole and neutron star systems and dwarf nova oscillations in cataclysmic variables.

Using an analytic, axisymmetric approach that includes general relativity, coupled to a condition for pair production deduced from simulations, we derive general results about the spatial distribution of pair-producing field lines over the pulsar polar cap. In particular, we show that pair production on magnetic field lines operates over only a fraction of the polar cap for an aligned rotator for general magnetic field configurations, assuming the magnetic field varies spatially on a scale that is larger than the size of the polar cap. We compare our result to force-free simulations of a pulsar with a dipole surface field and find excellent agreement. Our work has implications for first-principles simulations of pulsar magnetospheres, and for explaining observations of pulsed radio and high-energy emission.

We present PICsar -- a new Particle in Cell code geared towards efficiently simulating the magnetosphere of the aligned rotator. PICsar is a special relativistic, electromagnetic, charge conservative code that can be used to simulate arbitrary electromagnetics problems in axisymmetry. It features stretchable body-fitted coordinates that follow the surface of a sphere, simplifying the application of boundary conditions in the case of the aligned rotator; a radiation absorbing outer boundary, which allows a steady state to be set up dynamically and maintained indefinitely from transient initial conditions; and algorithms for injection of charged particles into the simulation domain. The code is parallelized using MPI and scales well to a large number of processors. We discuss the numerical methods used in PICsar and present tests of the code. In particular, we show that PICsar can accurately and efficiently simulate the magnetosphere of the aligned monopole rotator in the force free limit. We present simulations of the aligned dipole rotator in a forthcoming paper.

We study low frequency waves that propagate in a region of layered semi-convection. Layered semi-convection is predicted to be present in stellar and planetary interiors and can significantly modify the rate of thermal and compositional mixing. We derive a series of analytical dispersion relations for plane-parallel layered semi-convection in the Boussinesq approximation using a matrix transfer formalism. We find that like a continuously stratified medium, a semi-convective staircase -- in which small convective regions are separated by sharp density jumps -- supports internal gravity waves (g-modes). When the wavelength is much longer than the distance between semi-convective steps, these behave nearly like g-modes in a continuously stratified medium. However, the g-mode period spacing in a semi-convective region is systematically {\em smaller} than in a continuously stratified medium, and it decreases with decreasing mode frequency. When the g-mode wavelength becomes comparable to the distance between semi-convective steps, the g-mode frequencies deviate significantly from those of a continuously stratified medium (the frequencies are higher). G-modes with vertical wavelengths smaller than the distance between semi-convective steps are evanescent and do not propagate in the staircase. Thus, there is a lower cutoff frequency for a given horizontal wavenumber. We generalize our results to gravito-inertial waves relevant for rapidly rotating stars and planets. Finally, we assess the prospects for detecting layered semi-convection using astero/planetary seismology.

We perform 2.5D axisymmetric simulations of the pulsar magnetosphere (aligned dipole rotator) using the charge conservative, relativistic, electromagnetic particle in cell code PICsar. Particle in cell codes are a powerful tool to use for studying the pulsar magnetosphere, because they can handle the force-free and vacuum limits and provide a self-consistent treatment of magnetic reconnection. In the limit of dense plasma throughout the magnetosphere, our solutions are everywhere in the force-free regime except for dissipative regions at the polar caps, in the current layers, and at the Y-point. These dissipative regions arise self-consistently, since we do not have any explicit dissipation in the code. A minimum of $ \approx 15-20\%$ of the electromagnetic spindown luminosity is transferred to the particles inside 5 light cylinder radii. However, the particles can carry as much as $\gtrsim 50 \%$ of the spindown luminosity if there is insufficient plasma in the outer magnetosphere to screen the component of electric field parallel to the magnetic field. In reality, the component of the spindown luminosity carried by the particles could be radiated as gamma rays, but high-frequency synchrotron emission would need to be implemented as a sub-grid process in our simulations and is not present for the current suite of runs. The value of the spindown luminosity in our simulations is within $\approx 10\%$ of the force-free value, and the structure of the electromagnetic fields in the magnetosphere is on the whole consistent with the force-free model.

We perform global unstratified 3D magnetohydrodynamic simulations of an astrophysical boundary layer (BL) -- an interface region between an accretion disk and a weakly magnetized accreting object such as a white dwarf -- with the goal of understanding the effects of magnetic field on the BL. We use cylindrical coordinates with an isothermal equation of state and investigate a number of initial field geometries including toroidal, vertical, and vertical with zero net flux. Our initial setup consists of a Keplerian disk attached to a non-rotating star. In a previous work, we found that in hydrodynamical simulations, sound waves excited by shear in the BL were able to efficiently transport angular momentum and drive mass accretion onto the star. Here we confirm that in MHD simulations, waves serve as an efficient means of angular momentum transport in the vicinity of the BL, despite the magnetorotational instability (MRI) operating in the disk. In particular, the angular momentum current due to waves is at times larger than the angular momentum current due to MRI. Our results suggest that angular momentum transport in the BL and its vicinity is a global phenomenon occurring through dissipation of waves and shocks. This point of view is quite different from the standard picture of transport by a local anomalous turbulent viscosity. In addition to angular momentum transport, we also study magnetic field amplification within the BL. We find that the field is indeed amplified in the BL, but only by a factor of a few and remains subthermal.

The nature of angular momentum transport in the boundary layers of accretion disks has been one of the central and long-standing issues of accretion disk theory. In this work we demonstrate that acoustic waves excited by supersonic shear in the boundary layer serve as an efficient mechanism of mass, momentum and energy transport at the interface between the disk and the accreting object. We develop the theory of angular momentum transport by acoustic modes in the boundary layer, and support our findings with 3D hydrodynamical simulations, using an isothermal equation of state. Our first major result is the identification of three types of global modes in the boundary layer. We derive dispersion relations for each of these modes that accurately capture the pattern speeds observed in simulations to within a few percent. Second, we show that angular momentum transport in the boundary layer is intrinsically nonlocal, and is driven by radiation of angular momentum away from the boundary layer into both the star and the disk. The picture of angular momentum transport in the boundary layer by waves that can travel large distances before dissipating and redistributing angular momentum and energy to the disk and star is incompatible with the conventional notion of local transport by turbulent stresses. Our results have important implications for semianalytical models that describe the spectral emission from boundary layers.

We propose a new method to conserve the total energy to round-off error in grid-based codes for hydrodynamic simulations with self-gravity. A formula for the energy flux due to the work done by the the self-gravitational force is given, so the change in total energy can be written in conservative form. Numerical experiments with the code Athena show that the total energy is indeed conserved with our new algorithm and the new algorithm is second order accurate. We have performed a set of tests that show the numerical errors in the traditional, non-conservative algorithm can affect the dynamics of the system. The new algorithm only requires one extra solution of the Poisson equation, as compared to the traditional algorithm which includes self-gravity as a source term. If the Poisson solver takes a negligible fraction of the total simulation time, such as when FFTs are used, the new algorithm is almost as efficient as the original method. This new algorithm is useful in Eulerian hydrodynamic simulations with self-gravity, especially when results are sensitive to small energy errors, as for radiation pressure dominated flow.

Disk accretion onto a weakly magnetized central object, e.g. a star, is inevitably accompanied by the formation of a boundary layer near the surface, in which matter slows down from the highly supersonic orbital velocity of the disk to the rotational velocity of the star. We perform high resolution 2D hydrodynamical simulations in the equatorial plane of an astrophysical boundary layer with the goal of exploring the dynamics of non-axisymmetric structures that form there. We generically find that the supersonic shear in the boundary layer excites non-axisymmetric quasi-stationary acoustic modes that are trapped between the surface of the star and a Lindblad resonance in the disk. These modes rotate in a prograde fashion, are stable for hundreds of orbital periods, and have a pattern speed that is less than and of order the rotational velocity at the inner edge of the disk. The origin of these intrinsically global modes is intimately related to the operation of a corotation amplifier in the system. Dissipation of acoustic modes in weak shocks provides a universal mechanism for angular momentum and mass transport even in purely hydrodynamic (i.e. non-magnetized) boundary layers. We discuss the possible implications of these trapped modes for explaining the variability seen in accreting compact objects.

Disk accretion onto weakly magnetized astrophysical objects often proceeds via a boundary layer that forms near the object's surface, in which the rotation speed of the accreted gas changes rapidly. Here we study the initial stages of formation for such a boundary layer around a white dwarf or a young star by examining the hydrodynamical shear instabilities that may initiate mixing and momentum transport between the two fluids of different densities moving supersonically with respect to each other. We find that an initially laminar boundary layer is unstable to two different kinds of instabilities. One is an instability of a supersonic vortex sheet (implying a discontinuous initial profile of the angular speed of the gas) in the presence of gravity, which we find to have a growth rate of order (but less than) the orbital frequency. The other is a sonic instability of a finite width, supersonic shear layer, which is similar to the Papaloizou-Pringle instability. It has a growth rate proportional to the shear inside the transition layer, which is of order the orbital frequency times the ratio of stellar radius to the boundary layer thickness. For a boundary layer that is thin compared to the radius of the star, the shear rate is much larger than the orbital frequency. Thus, we conclude that sonic instabilities play a dominant role in the initial stages of nonmagnetic boundary layer formation and give rise to very fast mixing between disk gas and stellar fluid in the supersonic regime.

Collisions resulting in fragmentation are important in shaping the mass spectrum of minor bodies in the asteroid belt, the Kuiper belt, and debris disks. Models of fragmentation cascades typically find that in steady-state, the solution for the particle mass distribution is a power law in the mass. However, previous studies have typically assumed that the mass of the largest fragment produced in a collision with just enough energy to shatter the target and disperse half its mass to infinity is directly proportional to the target mass. We show that if this assumption is not satisfied, then the power law solution for the steady-state particle mass distribution is modified by a multiplicative factor, which is a slowly varying function of the mass. We derive analytic solutions for this correction factor and confirm our results numerically. We find that this correction factor proves important when extrapolating over many orders of magnitude in mass, such as when inferring the number of large objects in a system based on infrared observations. In the course of our work, we have also discovered an unrelated type of non-power law behavior: waves can persist in the mass distribution of objects even in the absence of upper or lower cutoffs to the mass distribution or breaks in the strength law.

We study the dynamics of large dust grains >1 micron with orbits outside of the heliosphere (beyond 250 AU). Motion of the Solar System through the interstellar medium (ISM) at a velocity of 26 km/s subjects these particles to gas and Coulomb drag (grains are expected to be photoelectrically charged) as well as the Lorentz force and the electric force caused by the induction electric field. We show that to zeroth order the combined effect of these forces can be well described in the framework of the classical Stark problem: particle motion in a Keplerian potential subject to an additional constant force. Based on this analogy, we elucidate the circumstances in which the motion becomes unbound, and show that under local ISM conditions dust grains smaller than ~100 microns originating in the Oort Cloud (e.g. in collisions of comets) beyond 10000 AU are ejected from the Solar System under the action of the electric force. Orbital motion of larger, bound grains is described analytically using the orbit-averaged Hamiltonian approach and consists of orbital plane precession at a fixed semi-major axis, accompanied by the periodic variations of the inclination and eccentricity (the latter may approach unity in some cases). A more detailed analysis of the combined effect of gas and Coulomb drag shows it is possible to reduce particle semi-major axes, but that the degree of orbital decay is limited (a factor of several at best) by passages through atomic and molecular clouds, which easily eject small particles.