Intermediate and massive stars drive fast and powerful isotropic winds that interact with the winds of nearby stars in star clusters and the surrounding interstellar medium (ISM). Wind-ISM collisions generate astrospheres around these stars that contain hot $T\sim 10^7$ K gas that adiabatically expands. As individual bubbles expand and collide they become unstable, potentially driving turbulence in star clusters. In this paper we use hydrodynamic simulations to model a densely populated young star cluster within a homogeneous cloud to study stellar wind collisions with the surrounding ISM. We model a mass-segregated cluster of 20 B-type young main sequence stars with masses ranging from 3--17 $M_{\odot}$. We evolve the winds for $\sim$11 kyrs and show that wind-ISM collisions and over-lapping wind-blown bubbles around B-stars mixes the hot gas and ISM material generating Kolmogorov-like turbulence on small scales early in its evolution. We discuss how turbulence driven by stellar winds may impact the subsequent generation of star formation in the cluster

Stellar feedback in the form of radiation pressure and magnetically-driven collimated outflows may limit the maximum mass that a star can achieve and affect the star-formation efficiency of massive pre-stellar cores. Here we present a series of 3D adaptive mesh refinement radiation-magnetohydrodynamic simulations of the collapse of initially turbulent, massive pre-stellar cores. Our simulations include radiative feedback from both the direct stellar and dust-reprocessed radiation fields, and collimated outflow feedback from the accreting stars. We find that protostellar outflows punches holes in the dusty circumstellar gas along the star's polar directions, thereby increasing the size of optically thin regions through which radiation can escape. Precession of the outflows as the star's spin axis changes due to the turbulent accretion flow further broadens the outflow, and causes more material to be entrained. Additionally, the presence of magnetic fields in the entrained material leads to broader entrained outflows that escape the core. We compare the injected and entrained outflow properties and find that the entrained outflow mass is a factor of $\sim$3 larger than the injected mass and the momentum and energy contained in the entrained material are $\sim$25% and $\sim$5% of the injected momentum and energy, respectively. As a result, we find that, when one includes both outflows and radiation pressure, the former are a much more effective and important feedback mechanism, even for massive stars with significant radiative outputs.

Star formation is a multi-scale, multi-physics problem ranging from the size scale of molecular clouds ($\sim$10s pc) down to the size scales of dense prestellar cores ($\sim$0.1 pc) that are the birth sites of stars. Several physical processes like turbulence, magnetic fields and stellar feedback, such as radiation pressure and outflows, are more or less important for different stellar masses and size scales. During the last decade a variety of technological and computing advances have transformed our understanding of star formation through the use of multi-wavelength observations, large scale observational surveys, and multi-physics multi-dimensional numerical simulations. Additionally, the use of synthetic observations of simulations have provided a useful tool to interpret observational data and evaluate the importance of various physical processes on different scales in star formation. Here, we review these recent advancements in both high- ($M \gtrsim 8 \, M_{\rm \odot}$) and low-mass star formation.

Recent observations have suggested that circumstellar disks may commonly form around young stellar objects. Although the formation of circumstellar disks can be a natural result of the conservation of angular momentum in the parent cloud, theoretical studies instead show disk formation to be difficult from dense molecular cores magnetized to a realistic level, owing to efficient magnetic braking that transports a large fraction of the angular momentum away from the circumstellar region. We review recent progress in the formation and early evolution of disks around young stellar objects of both low-mass and high-mass, with an emphasis on mechanisms that may bridge the gap between observation and theory, including non-ideal MHD effects and asymmetric perturbations in the collapsing core (e.g., magnetic field misalignment and turbulence). We also address the associated processes of outflow launching and the formation of multiple systems, and discuss possible implications in properties of protoplanetary disks.

Using numerical hydrodynamics calculations and a novel method for densely sampling parameter space, we measure the accretion and torque on a binary system from a circumbinary disk. In agreement with previous studies, we find that the net torque on the binary is positive for mass ratios close to unity, and that accretion always drives the binary towards equal mass. Accretion variability depends sensitively on the numerical sink prescription, but the torque and relative accretion onto each component do not depend on the sink timescale. Positive torque and highly variable accretion occurs only for mass ratios greater than around $0.05$. This means that for mass ratios below $0.05$, the binary would migrate inward until the secondary accreted sufficient mass, after which it would execute a U-turn and migrate outward. We explore a range of viscosities, from $\alpha = 0.03$ to $\alpha = 0.15$, and find that this outward torque is proportional to the viscous torque, simply proportional to viscosity in this range. Dependence of accretion and torque on mass ratio is explored in detail, densely sampling mass ratios between $0.01$ and unity. For mass ratio $q > 0.6$, accretion variability is found to exhibit a distinct sawtooth pattern, typically with a five-orbit cycle that provides a "smoking gun" prediction for variable quasars observed over long periods, as a potential means to confirm the presence of a binary.

Similar to their low-mass counterparts, massive stars likely form via the collapse of pre-stellar molecular cores. Recent observations suggest that most massive cores are subvirial (i.e., not supported by turbulence) and therefore are likely unstable to gravitational collapse. Here we perform radiation hydrodynamic simulations to follow the collapse of turbulent massive pre-stellar cores with subvirial and virialized initial conditions to explore how their dynamic state affects the formation of massive stars and core fragmentation into companion stars. We find that subvirial cores undergo rapid monolithic collapse resulting in higher accretion rates at early times as compared to the collapse of virialized cores that have the same physical properties. In contrast, we find that virialized cores undergo a slower, gradual collapse and significant turbulent fragmentation at early times resulting in numerous companion stars. In the absence of strong magnetic fields and protostellar outflows we find that the faster growth rate of massive stars that are born out of subvirial cores leads to an increase in the radiative heating of the core thereby further suppressing fragmentation at early times when turbulent fragmentation occurs for virialized cores. Regardless of initial condition, we find that the massive accretion disks that form around massive stars dominant the accretion flow onto the star at late times and eventually become gravitationally unstable and fragment to form companion stars at late times.

Software that processes real-world data or that models a physical system must have some way of managing units. While simple approaches like the understood convention that all data are in a unit system (such as the MKS SI unit system) do work in practice, they are fraught with possible sources of error both by developers and users of the software. In this paper we present unyt, a Python library based on NumPy and SymPy for handling data that has units. It is designed both to aid quick interactive calculations and to be tightly integrated into a larger Python application or library. We compare unyt with two other Python libraries for handling units, Pint and astropy.units, and find that unyt is faster, has higher test coverage, and has fewer lines of code.

The pressure exerted by massive stars' radiation fields is an important mechanism regulating their formation. Detailed simulation of massive star formation therefore requires an accurate treatment of radiation. However, all published simulations have either used a diffusion approximation of limited validity; have only been able to simulate a single star fixed in space, thereby suppressing potentially-important instabilities; or did not provide adequate resolution at locations where instabilities may develop. To remedy this we have developed a new, highly accurate radiation algorithm that properly treats the absorption of the direct radiation field from stars and the re-emission and processing by interstellar dust. We use our new tool to perform three-dimensional radiation-hydrodynamic simulations of the collapse of massive pre-stellar cores with laminar and turbulent initial conditions and properly resolve regions where we expect instabilities to grow. We find that mass is channeled to the stellar system via gravitational and Rayleigh-Taylor (RT) instabilities, in agreement with previous results using stars capable of moving, but in disagreement with methods where the star is held fixed or with simulations that do not adequately resolve the development of RT instabilities. For laminar initial conditions, proper treatment of the direct radiation field produces later onset of instability, but does not suppress it entirely provided the edges of radiation-dominated bubbles are adequately resolved. Instabilities arise immediately for turbulent pre-stellar cores because the initial turbulence seeds the instabilities. Our results suggest that RT features are significant and should be present around accreting massive stars throughout their formation.

We present a highly-parallel multi-frequency hybrid radiation hydrodynamics algorithm that combines a spatially-adaptive long characteristics method for the radiation field from point sources with a moment method that handles the diffuse radiation field produced by a volume-filling fluid. Our Hybrid Adaptive Ray-Moment Method (HARM$^2$) operates on patch-based adaptive grids, is compatible with asynchronous time stepping, and works with any moment method. In comparison to previous long characteristics methods, we have greatly improved the parallel performance of the adaptive long-characteristics method by developing a new completely asynchronous and non-blocking communication algorithm. As a result of this improvement, our implementation achieves near-perfect scaling up to $\mathcal{O}(10^3)$ processors on distributed memory machines. We present a series of tests to demonstrate the accuracy and performance of the method.

Star clusters larger than $\sim 10^{3}$ $M_\odot$ contain multiple hot stars that launch fast stellar winds. The integrated kinetic energy carried by these winds is comparable to that delivered by supernova explosions, suggesting that at early times winds could be an important form of feedback on the surrounding cold material from which the star cluster formed. However, the interaction of these winds with the surrounding clumpy, turbulent, cold gas is complex and poorly understood. Here we investigate this problem via an accounting exercise: we use empirically determined properties of four well-studied massive star clusters to determine where the energy injected by stellar winds ultimately ends up. We consider a range of kinetic energy loss channels, including radiative cooling, mechanical work on the cold interstellar medium, thermal conduction, heating of dust via collisions by the hot gas, and bulk advection of thermal energy by the hot gas. We show that, for at least some of the clusters, none of these channels can account for more than a small fraction of the injected energy. We suggest that turbulent mixing at the hot-cold interface or physical leakage of the hot gas from the HII region can efficiently remove the kinetic energy injected by the massive stars in young star clusters. Even for the clusters where we are able to account for all the injected kinetic energy, we show that our accounting sets strong constraints on the importance of stellar winds as a mechanism for feedback on the cold interstellar medium.

Recent studies have claimed the existence of very massive stars (VMS) up to 300 solar masses in the local Universe. As this finding may represent a paradigm shift for the canonical stellar upper-mass limit of 150 Msun, it is timely to discuss the status of the data, as well as the far-reaching implications of such objects. We held a Joint Discussion at the General Assembly in Beijing to discuss (i) the determination of the current masses of the most massive stars, (ii) the formation of VMS, (iii) their mass loss, and (iv) their evolution and final fate. The prime aim was to reach broad consensus between observers and theorists on how to identify and quantify the dominant physical processes.

The physical mechanisms that set the initial rotation rates in massive stars are a crucial unknown in current star formation theory. Observations of young, massive stars provide evidence that they form in a similar fashion to their low-mass counterparts. The magnetic coupling between a star and its accretion disk may be sufficient to spin down low-mass pre-main sequence (PMS) stars to well below breakup at the end stage of their formation when the accretion rate is low. However, we show that these magnetic torques are insufficient to spin down massive PMS stars due to their short formation times and high accretion rates. We develop a model for the angular momentum evolution of stars over a wide range in mass, considering both magnetic and gravitational torques. We find that magnetic torques are unable to spin down either low or high mass stars during the main accretion phase, and that massive stars cannot be spun down significantly by magnetic torques during the end stage of their formation either. Spin-down occurs only if massive stars' disk lifetimes are substantially longer or their magnetic fields are much stronger than current observations suggest.

We investigate potential models that could explain why multiple proto-stellar systems predominantly show single jets. During their formation, stars most frequently produce energetic outflows and jets. However, binary jets have only been observed in a very small number of systems. We model numerically 3D binary jets for various outflow parameters. We also model the propagation of jets from a specific source, namely L1551 IRS 5, known to have two jets, using recent observations as constraints for simulations with a new MHD code. We examine their morphology and dynamics, and produce synthetic emission maps. We find that the two jets interfere up to the stage where one of them is almost destroyed or engulfed into the second one. We are able to reproduce some of the observational features of L1551 such as the bending of the secondary jet. While the effects of orbital motion are negligible over the jets dynamical timeline, their interaction has significant impact on their morphology. If the jets are not strictly parallel, as in most observed cases, we show that the magnetic field can help the collimation and refocusing of both of the two jets.

NGC 7538 IRS1 is a high-mass (approx. 30 M_sun) protostar with a CO outflow, an associated UCHII region, and a linear methanol maser structure, which might trace a Keplerian-rotating circumstellar disk. The directions of the various associated axes are misaligned with each other. We investigate the near-infrared morphology of the source to clarify the relations among the various axes. K'-band bispectrum speckle interferometry was performed at two 6-meter-class telescopes -- the BTA 6m telescope and the 6.5m MMT. Complementary IRAC images from the Spitzer Space Telescope Archive were used to relate the structures detected with the outflow at larger scales. High-dynamic range images show fan-shaped outflow structure in which we detect 18 stars and several blobs of diffuse emission. We interpret the misalignment of various outflow axes in the context of a disk precession model, including numerical hydrodynamic simulations of the molecular emission. The precession period is approx. 280 years and its half-opening angle is 40 degrees. A possible triggering mechanism is non-coplanar tidal interaction of an (undiscovered) close companion with the circumbinary protostellar disk. Our observations resolve the nearby massive protostar NGC 7538 IRS2 as a close binary with separation of 195 mas. We find indications for shock interaction between the outflow activities in IRS1 and IRS2. Indications of outflow precession have been discovered to date in a number of massive protostars, all with large precession angles 20--45 degrees. This might explain the difference between the outflow widths in low- and high-mass stars and add support to a common collimation mechanism.

We review the known properties of molecular outflows from low- and high-mass young stars. General trends among outflows are identified, and the most recent studies on the morphology, kinematics, energetics, and evolution of molecular outflows are discussed, focusing on results from high-resolution millimeter observations. We review the existing four broad classes of outflow models and compare numerical simulations with the observational data. A single class of models cannot explain the range of morphological and kinematic properties that are observed, and we propose a possible solution. The impact of outflows on their cloud is examined, and we review how outflows can disrupt their surrounding environment, through the clearing of gas and the injection of momentum and energy onto the gas at distances from their powering sources from about 0.01 to a few pc. We also discuss the effects of shock-induced chemical processes on the ambient medium, and how these processes may act as a chemical clock to date outflows. Lastly, future outflow research with existing and planned millimeter and submillimeter instruments is presented.

We present hydrodynamic simulations of molecular outflows driven by jets with a long period of precession, motivated by observations of arc-like features and S-symmetry in outflows associated with young stars. We simulate images of not only H2 vibrational and CO rotational emission lines, but also of atomic emission. The density cross section displays a jaw-like cavity, independent of precession rate. In molecular hydrogen, however, we find ordered chains of bow shocks and meandering streamers which contrast with the chaotic structure produced by jets in rapid precession. A feature particularly dominant in atomic emission is a stagnant point in the flow that remains near the inlet and alters shape and brightness as the jet skims by. Under the present conditions, slow jet precession yields a relatively high fraction of mass accelerated to high speeds, as also attested to in simulated CO line profiles. Many outflow structures, characterised by HH222 (continuous ribbon), HH240 (asymmetric chains of bow shocks) and RNO43N (protruding cavities), are probably related to the slow precession model.

The kinematic and radiative power of molecular jets is expected to change as a protostar undergoes permanent or episodal changes in the rate at which it accretes. We study here the consequences of evolving jet power on the spatial and velocity structure, as well as the fluxes, of molecular emission from the bipolar outflow. We consider a jet of rapidly increasing density and a jet in which the mass input is abruptly cut off. We perform three dimensional hydrodynamic simulations with atomic and molecular cooling and chemistry. In this work, highly collimated and sheared jets are assumed. We find that position-velocity diagrams, velocity-channel maps and the relative H$_2$ and CO fluxes are potentially the best indicators of the evolutionary stage. In particular, the velocity width of the CO lines may prove most reliable although the often-quoted mass-velocity power-law index is probably not. We demonstrate how the relative H$_2$ 1--0 S(1) and CO J=1--0 fluxes evolve and apply this to interpret the phase of several outflows.

We present the results from three dimensional hydrodynamical simulations of decaying high-speed turbulence in dense molecular clouds. We compare our results, which include a detailed cooling function, molecular hydrogen chemistry and a limited C and O chemistry, to those previously obtained for decaying isothermal turbulence. After an initial phase of shock formation, power-law decay regimes are uncovered, as in the isothermal case. We find that the turbulence decays faster than in the isothermal case because the average Mach number remains higher, due to the radiative cooling. The total thermal energy, initially raised by the introduction of turbulence, decays only a little slower than the kinetic energy. We discover that molecule reformation, as the fast turbulence decays, is several times faster than that predicted for a non-turbulent medium. This is caused by moderate speed shocks which sweep through a large fraction of the volume, compressing the gas and dust. Through reformation, the molecular density and molecular column appear as complex patterns of filaments, clumps and some diffuse structure. In contrast, the molecular fraction has a wider distribution of highly distorted clumps and copious diffuse structure, so that density and molecular density are almost identically distributed during the reformation phase. We conclude that molecules form in swept-up clumps but effectively mix throughout via subsequent expansions and compressions.

Two three-dimensional magnetohydrodynamical simulations of strongly magnetized conical jets, one with a poloidal and one with a helical magnetic field, have been performed. In the poloidal simulation a significant sheath (wind) of magnetized moving material developed and partially stabilized the jet to helical twisting. The fundamental pinch mode was not similarly affected and emission knots developed in the poloidal simulation. Thus, astrophysical jets surrounded by outflowing winds could develop knotty structures along a straight jet triggered by pinching. Where helical twisting dominated the dynamics, magnetic field orientation along the line-of-sight could be organized by the toroidal flow field accompanying helical twisting. On astrophysical jets such structure could lead to a reversal of the direction of Faraday rotation in adjacent zones along a jet. Theoretical analysis showed that the different dynamical behavior of the two simulations could be entirely understood as a result of dependence on the velocity shear between jet and wind which must exceed a surface Alfven speed before the jet becomes unstable to helical and higher order modes of jet distortion.

Three dimensional numerical simulations of the response of a Lorentz factor 2.5 relativistic jet to precession at three different frequencies have been performed. Low, moderate and high precession frequencies have been chosen relative to the maximally unstable frequency predicted by a Kelvin-Helmholtz stability analysis. Transverse motion and velocity decreases as the precession frequency increases. Although helical displacement of the jet decreases in amplitude as the precession frequency increases, a helical shock is generated in the medium external to the jet at all precession frequencies. Complex pressure and velocity structure inside the jet is shown to be produced by a combination of the helical surface and first body modes predicted by a normal mode analysis of the relativistic hydrodynamic equations. The surface and first body mode have different wave speed and wavelength, are launched in phase by the periodic precession, and exhibit beat patterns in synthetic emission images. Wave (pattern) speeds range from $0.41c$ to $0.86c$ but beat patterns remain stationary. Thus, we find a mechanism that can produce differentially moving and stationary features in the jet.

We extend investigations of mass entrainment by jets, which previously have focused on cylindrical supermagnetosonic jets and expanding trans-Alfv\'enic jets, to a set of expanding supermagnetosonic jets. We precess these jets at the origin to excite the helical mode of the Kelvin-Helmholtz (or KH) instability, in order to compare the results with predictions from linear stability analysis. We analyze this simulation set for the spatial development of magnetized mass, which we interpret as jet plus entrained, initially unmagnetized external mass. As with the previous simulation sets, we find that the growth of magnetized mass is associated with the growth of the KH instability through linear, nonlinear, and saturated stages and with the expansion of magnetized material in simulated observations of the jet. From comparison of measured wavelengths and wave speeds with the predictions from linear stability analysis, we see evidence that the KH instability is the primary cause for mass entrainment in these simulations, and that the expansion reduces the rate of mass entrainment. This reduced rate can be observed as a somewhat greater distance between the two transition points separating the three stages of expansion.

We investigate the growth of jet plus entrained mass in simulations of supermagnetosonic cylindrical and expanding jets. The entrained mass spatially grows in three stages: from an initially slow spatial rate to a faster rate and finally at a flatter rate. These stages roughly coincide with the similar rates of expansion in simulated radio intensity maps, and also appear related to the growth of the Kelvin-Helmholtz instability through linear, nonlinear, and saturated regimes. In the supermagnetosonic cylindrical jets, we found that a jet with an embedded primarily toroidal magnetic field is more stable than a jet with a primarily axial magnetic field. Also, pressure-matched expanding jets are more stable and entrain less mass than cylindrical jets with equivalent inlet conditions.

The amplitude of jet distortions and accompanying pressure and velocity fluctuations resulting from Kelvin-Helmholtz instability of three dimensional relativistic jets are explored. The effect of instability on jets as they accelerate from sub- to super-Alfv\'enic speeds is explored and a relativistic stabilization mechanism for trans-Alfv\'enic jets is proposed. The level to which asymmetric instabilities on supermagnetosonic relativistic jets will grow is predicted theoretically and a Doppler boosted ``apparent'' emissivity is computed. Effects due to helically twisted filamentary structure produced by asymmetric modes of instability should be readily observable on relativistic jets.

Three-dimensional magnetohydrodynamical simulations of strongly magnetized ``light'' conical jets have been performed. An investigation of the transition from sub-Alfv\'enic to super-Alfv\'enic flow has been made for nearly poloidal and for helical magnetic fields. The jets are stable to asymmetric modes of jet distortion provided they are sub-Alfv\'enic over most of their interior but destabilize rapidly when they become on average super-Alfv\'enic. The jets are precessed at the origin and the resulting small amplitude azimuthal motion is communicated down the jet to the Alfv\'en point where it couples to a slowly moving and rapidly growing helical twist. Significant jet rotation can contribute to destabilization via increase in the velocity shear between the jet and the external medium. Destabilization is accompanied by significant mass entrainment and the jets slow down significantly as denser external material is entrained. Synchrotron intensity images satisfactorily reveal large scale helical structures but have trouble distinguishing a large amplitude elliptical jet distortion that appears as an apparent pinching in an intensity image. Smaller scale jet distortions are not clearly revealed in intensity images, largely as a result of the relatively small total pressure variations that accompany destabilization and growing distortions. Fractional polarization is high as a result of the strong ordered magnetic fields except where the intensity image suggests cancellation of polarization vectors by integration through twisted structures.

We compare results from a relativistic and a nonrelativistic set of 2D axisymmetric jet simulations. For a set of five relativistic simulations that either increase the Lorentz factor or decrease the adiabatic index we compute nonrelativistic simulations with equal useful power or thrust. We examine these simulations for morphological and dynamical differences, focusing on the velocity field, the width of the cocoon, the age of the jets, and the internal structure of the jet itself. The primary result of these comparisons is that the velocity field of nonrelativistic jet simulations cannot be scaled up to give the spatial distribution of Lorentz factors seen in relativistic simulations. Since the local Lorentz factor plays a major role in determining the total intensity for parsec scale extragalactic jets, this suggests that a nonrelativistic simulation cannot yield the proper intensity distribution for a relativistic jet. Another general result is that each relativistic jet and its nonrelativistic equivalents have similar ages (in dynamical time units, = R/a_a, where R is the initial radius of a cylindrical jet and a_a is the sound speed in the ambient medium). In addition to these comparisons, we have completed four new relativistic simulations to investigate the effect of varying thermal pressure on relativistic jets. The simulations generally confirm that faster (larger Lorentz factor) and colder jets are more stable, with smaller amplitude and longer wavelength internal variations. The apparent stability of these jets does not follow from linear normal mode analysis, which suggests that there are available growing Kelvin-Helmholtz modes. (Abridged.)

A set of 3D MHD simulations of magnetized jets has been performed. The jets contain an equipartition primarily poloidal magnetic field and the effect of jet density on jet dynamics and structure is evaluated. The jet is precessed at the origin to excite Kelvin-Helmholtz unstable helical modes. We extensively compare the structure in these simulations with linear stability theory. The jet that is dense with respect to the external medium develops a high speed core surrounded by a less dense sheath consisting of slower moving jet fluid. These simulations suggest that extended extragalactic jets propagate to such large distances because they are surrounded by a lobe or cocoon whose density is less than the jet density. (Abridged abstract.)

In a previous paper, we calculated numerical hydrodynamic models of the interstellar medium in the Galaxy, which suggested that hot gas (T $\ge$ 3 $\times$ 10$^5$ K) has a filling factor near 50\% in the midplane (Rosen \& Bregman 1995, ApJ 440, 634). Here we extend the work to examine the X-ray emission characteristics of the best model and calculate a variety of observable measures for comparison with the observed soft X-ray background. We calculate the X-ray intensities, spectra, and hardness ratios in two soft X-ray bands as Galactic latitude and cool gas column are varied, and make comparisons with observational data. The observed distribution of N$_{H I}$ is similar to that seen from one of the simulated bubbles, with the important exception that the minimum hydrogen column in the model is too low. Also, there are clear examples of anticorrelations between H I and X-ray emission as well as correlations between H I and X-rays, which is caused by an increased emission measure as a shock enters a cool gas region. Statistically, anticorrelations are slightly more common than correlations. X-ray spectra are calculated from the models and these reveal that for observations to have strong diagnostic power in probing the hot ISM, a spectral resolution of E/$\Delta$E $>$ 30 is required. This analysis demonstrates several shortcomings of our model, which may be solved by future simulations that include magnetic fields.