© The Authors 2024

1 Introduction

Star formation has been studied over decades both theoretically and numerically on large spatial scales, the key motivation being to determine how global redistribution of angular momentum results in the formation of stars, disks, or jets and set the initial conditions for planet formation (see Larson 1969; Penston 1969 and e.g., Pattle et al. 2023; Pineda et al. 2023; Tsukamoto et al. 2023 for recent reviews). The field has been stimulated by instruments capable of resolving the largest scales of these young stellar systems in ever-finer detail in several wavelength domains, such as the Atacama Large Millimeter Array (ALMA), the Very Large Telescope (VLT) or the James Webb Space Telescope (JWST). The importance of these questions has perhaps overshadowed that of the physical processes taking place simultaneously on a small scale. One question, for example, is whether it is possible for dust to start concentrating as soon as the cloud collapses, thereby promoting highly efficient planet formation later on.

The local box formalism provides a simple treatment of the hydrodynamics of a spherically symmetric system undergoing contraction. The global effects of curvature and acceleration due to collapse are consistently incorporated into the local treatment of the problem by source terms that can be time-dependent, while allowing for the adoption of periodic boundaries as the system is spatially homogeneous. Numerical simulations using the local approximation approach play an important role in the study of local gas dynamics and instabilities in astrophysics. The best-known example is the rediscovery of the magnetorotational instability (MRI) with the local shear box approximation, which was first applied to magnetohydrodynamic (MHD) simulations by Hawley et al. (1995). Local models have been used for (isotropic) cosmological collapses (Robertson & Goldreich 2012). Anisotropy in collapse and expansion was taken into account in the study of starless molecular core collapse (Toci et al. 2018), solar winds (Velli et al. 1992; Grappin et al. 1993; Grappin & Velli 1996; Tenerani & Velli 2017; Shi et al. 2020; Squire et al. 2020) (among which Tenerani & Velli 2017 also consider nonuniform expansion), anisotropy-driven plasma kinetic instabilities (Sironi & Narayan 2015; Bott et al. 2021), and the cosmic ray pressure anisotropy instability (Sun & Bai 2023).

Recently, Lynch & Laibe (2023) (hereafter LL23) developed a new local model for the spherical collapse or expansion problem, allowing for a nonuniform collapse or expansion. This model considers a radially varying flow, as in Tenerani & Velli (2017), but also allows for the time-dependence of the background flow. It also uses a different treatment of pressure gradient than the expanding box models, in an approach that is more similar to the distorted shearing box models of Ogilvie & Latter (2013); Ogilvie & Barker (2014). This avoids spurious instabilities from the inclusion of formally subdominant terms (Latter & Papaloizou 2017), allowing for the application of periodic or shear-periodic boundary conditions with a Cartesianlike geometry. In this study, we incorporate LL23’s model into ATHENA++, a high-order Godunov MHD code (Stone et al. 2020). Other local models have already been implemented in ATHENA++, including Squire et al. (2020) for the (uniformly) expanding box, and MHD-PIC (Sun & Bai 2023), also allowing for an anisotropic collapse or expansion. However, most of these previous numerical implementations focus on the study of MHD properties and magnetized waves, in the context of stellar winds or cosmic rays, and the hydrodynamic properties of the expanding or collapsing box have not been comprehensively tested so far. Here, we aim to incorporate the local model of LL23, with benchmark tests to thoroughly investigate and validate the hydrodynamic properties therein.

This article is organized as follows. We first describe the formulation of the local model for the collapsing box in Sec. 2. In Sect. 3, we describe the numerical implementation, especially the modification of the Riemann solver in order to incorporate the local model. We then conduct a series of test cases including nonlinear solutions and linear perturbations in Sec. 4, as well as numerical studies. We further discuss possible applications and perspectives, and give future perspectives in Sec. 5. We summarize and conclude in Sec. 6.

2 Equations of the local model

We started by considering an axisymmetric spherical gas cloud collapsing or expanding with a purely radial velocity, U_bg(r, t), that varies with both radius, r, and time, t. We then considered a local neighborhood of a reference point located at a radius, ℛ(t), in the spherical gas cloud, described by a local Cartesian coordinate system comoving with the background flow. The local coordinates are given by $$x=\varphi ,y=\theta -\frac{\pi }{2},z=\frac{r-\mathcal{R}(t)}{L_z(t)}\text{, }$$(1)

where ϕ and θ are the longitude and latitude on the reference sphere of radius, ℛ, and r − ℛ(t) is the local radial distance relative to the global, ℛ(t), of the reference frame. We note that x, y, z have dimensions of angles in this formalism, and that they denote two azimuthal and radial directions, respectvely. L_z is a characteristic vertical length scale that is set by the background flow profile, which can be approximated by $$U_{bg}(r,t)=U_0(t)+U_{R0}(t)(r-\mathcal{R}(t)).$$(2)

The Lagrangian time derivative in the local model is $$D={\mathrm{\partial }}_t+v_i{\mathrm{\partial }}_i,$$(3)

while the divergence of the velocity of the radial background flow is $$\mathrm{\Delta }=\frac{2U_0}{\mathcal{R}}+U_{R0}.$$(4)

The flow within the local box model represents local, nonlinear perturbations to the background flow. Following the derivation in Sect. 2 of LL23, in the local frame, the continuity, momentum, and energy equations are given in the stretched Cartesian coordinate system by $$D\rho =-\rho [\mathrm{\Delta }+{\mathrm{\partial }}_x{v}_x+{\mathrm{\partial }}_y{v}_y+{\mathrm{\partial }}_z{v}_z],$$(5) $$D{v}_x+\frac{2U_0}{\mathcal{R}}{v}_x=-\frac{1}{\rho {\mathcal{R}}^2}{\mathrm{\partial }}_{x}p,$$(6) $$D{v}_y+\frac{2U_0}{\mathcal{R}}{v}_y=-\frac{1}{\rho {\mathcal{R}}^2}{\mathrm{\partial }}_{y}p,$$(7) $$D{v}_z+2v_z{U}_{R0}=-\frac{1}{\rho L_z^2}{\mathrm{\partial }}_{z}p,$$(8) $$\begin{array}{rl}& {\mathrm{\partial }}_t\left(\rho {\mathcal{E}}_{\text{rell }}\right)+{\mathrm{\partial }}_i\left[v_i(\rho {\mathcal{E}}_{\text{rell }}+p)\right]\\ & =-\rho [U_R{L}_z^2{v}_z^2+\frac{U_0}{\mathcal{R}}{\mathcal{R}}^2(v_x^2+v_y^2)]-(p+\rho {\mathcal{E}}_{\text{rell }})\mathrm{\Delta }\end{array},$$(9)

where ρ is the density of the gas, υ_i is the relative velocity field in the local frame¹, p is the pressure of the gas, and $\mathcal{E}$_rell is the energy density of the relative motion in this system, (re-)defined as $${\mathcal{E}}_{\text{rell }}=\frac{1}{2}{L}_z^2{v}_z^2+\frac{1}{2}{\mathcal{R}}^2(t)(v_x^2+v_y^2)+\epsilon ,$$(10)

which is the sum of the kinetic energy density of relative motion in the local frame, ${\mathcal{E}}_{k,\text{ rell }}\equiv \frac{1}{2}{L}_z^2{v}_z^2+\frac{1}{2}{\mathcal{R}}^2(t)(v_x^2+v_y^2)$, and the specific internal energy, ϵ, determined by the equation of state (EoS). For an isothermal EoS, $p=\rho c_{\mathrm{s}}^2$, where c_s is the sound speed, and the energy Eq. (9) is dropped from the system. For an adiabatic EoS, $\epsilon =\frac{p}{\rho (\gamma -1)}$, where γ is the ratio of specific heats, assumed to be constant.

The right-hand side of Eqs. (6)–(8) can be interpreted as an anisotropic sound speed in the local coordinate. We define two rescaled sound speeds in the horizontal and vertical directions separately: c_s,H ≡ c_s/ℛ and c_s,V = c_s/L_z.

Equations (5)–(9) can be written in a compact form, $$\frac{\mathrm{\partial }\mathbf{q}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\mathbf{f}}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\mathbf{g}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\mathbf{h}}{\mathrm{\partial }z}=\mathbf{S}(x,y,z,t),$$(11)

where the vector, q, of conservative variables is defined according to $$\mathbf{q}=\left(\begin{array}{l}{q}_1\\ q_2\\ q_3\\ q_4\\ q_5\end{array}\right)\equiv \left(\begin{array}{c}\rho \\ \rho v_x\\ \rho v_y\\ \rho v_z\\ \rho {\mathcal{E}}_{\text{rell }}\end{array}\right),$$(12)

f, g, and h are the fluxes of the associated conservative variables in the x, y, and z directions, respectvely: $$\mathbf{f}\equiv \left(\begin{array}{c}\rho v_x\\ \rho v_x^2+p/{\mathcal{R}}^2\\ \rho v_x{v}_y\\ \rho v_x{v}_z\\ (\rho {\mathcal{E}}_{\text{rell }}+p)v_x\end{array}\right),\mathbf{g}\equiv \left(\begin{array}{c}\rho v_y\\ \rho v_y{v}_x\\ \rho v_y^2+p/{\mathcal{R}}^2\\ \rho v_y{v}_z\\ (\rho {\mathcal{E}}_{\text{rell }}+p)v_y\end{array}\right),$$(13) $$\mathbf{h}\equiv \left(\begin{array}{c}\rho v_z\\ \rho v_z{v}_x\\ \rho v_z{v}_y\\ \rho v_z^2+p/L_z^2\\ (\rho {\mathcal{E}}_{\text{rell }}+p)v_z\end{array}\right).$$(14)

The vector of source terms, denoted as S, is $$\mathbf{S}\equiv \left(\begin{array}{c}-\rho \mathrm{\Delta }\\ -2\rho v_x\frac{U_0}{{\mathcal{R}}_0}-\rho \mathrm{\Delta }{v}_x\\ -2\rho v_y\frac{U_0}{\mathcal{R}}-\rho \mathrm{\Delta }{v}_y\\ -2\rho v_z{U}_{R0}-\rho \mathrm{\Delta }{v}_z\\ -\rho [U_R{L}_z^2{v}_z^2+\frac{U_0}{\mathcal{R}}{\mathcal{R}}^2(v_x^2+v_y^2)]-(p+\rho {\mathcal{E}}_{\text{rell }})\mathrm{\Delta }\end{array}\right).$$(15)

Again, for an isothermal EoS, the fifth component of the vectors q, f, g, h, and S are dropped, as the energy equation is not solved for this case.

A more intuitive illustration of the local model geometry is presented in Fig. 1. As the collapse proceeds, the local domain is compressed in volume, leading to an increase in gas density. This is reflected in the source term in the continuity Eq. (5). The angular momentum conservation law leads to amplification of horizontal flow velocity in the local frame during the collapse, which is represented by the source terms in the momentum equations on the horizontal directions (Eqs. (6) and (7)). These features along with other characteristic behaviors of the local model are further illustrated and tested in Sec. 4.

Fig. 1 Cartoon illustration of the local collapsing box model. On the left, the global geometry depicts a domain moving radially due to background flow. The volume of the domain varies as it approaches the center of the spherical cloud, leading to changes in density. A compressing domain leads to a density increase, as is shown by the progressively darker blue color in the sequence. The horizontal domain size is proportional to ℛ(t), establishing the horizontal length scale in the local model (middle panel). The vertical length scale, Lz (t), is determined by the distance between the domain’s inner and outer radial boundaries, which is set by the background flow. Inside each box, a waveform signal illustrates how certain physical quantities evolve in the local model. The width of the waveform represents quantities that remains fixed in the local frame (e.g., wavelength). The amplitude of the waveform represents specific physical quantities that are expected to be constant in the global geometry throughout the collapse (e.g., angular momentum, isothermal physical sound speed). In the local reference frame (right), the corresponding quantities such as the shear flow velocity or the effective sound speed are amplified as the collapse proceeds and the length scale shrinks in the domain, while quantities such as the wavelength remain the same. This phenomenon is visually demonstrated by the vertically stretched waveforms in the sequence. For a more accurate and detailed understanding of this process, one should refer to the equations in Sec. 2.

3 Numerical implementation

We next translated the local model of Sec. 2 into a Godunov scheme, and implemented it in the numerical code ATHENA++. Finite volume discretization of the compact form of the hydro-dynamic equations Eq. (11) gives the time integration step as a set of ordinary differential equations, $$\begin{array}{rl}{\mathbf{q}}_{i,j,k}^{n+1}=& {\mathbf{q}}_{i,j,k}^n-\frac{\delta t}{\delta x}({\mathbf{f}}_{i+1/2,j,k}^{n+1/2}-{\mathbf{f}}_{i-1/2,j,k}^{n+1/2})\\ & -\frac{\delta t}{\delta y}({\mathbf{g}}_{i,j+1/2,k}^{n+1/2}-{\mathbf{g}}_{i,j-1/2,k}^{n+1/2})\\ & -\frac{\delta t}{\delta z}({\mathbf{h}}_{i,j,k+1/2}^{n+1/2}-{\mathbf{h}}_{i,j,k-1/2}^{n+1/2})\\ & +\delta t{\mathbf{S}}_{i,j,k}^n\\ \equiv & {\tilde{\mathbf{q}}}_{i,j,k}^n+\delta t{\mathbf{S}}_{i,j,k}^n,\end{array}$$(16)

where subscript indices (i, j, k) denote numerical cells centered at positions (x_i, y_j, z_k), superscripts, n, denotes the index for discrete time, δx, δy, and δ_z are discrete grid sizes, and δ_t is the time step. The vectors ${\mathbf{q}}_{i,j,k}^n$ and ${\mathbf{S}}_{i,j,k}^n$ correspond to volume-averaged q and S over the numerical cells. The vectors ${\mathbf{f}}_{i+1/2,j,k}^{n+1/2},{\mathbf{g}}_{i+1/2,j,k}^{n+1/2},$, and ${\mathbf{h}}_{i+1/2,j,k}^{n+1/2}$ correspond to time-averaged f, g, and h, respectvely. The half-integer subscripts on the flux vectors refer to the edges of computational cells, while the half-integer superscripts denote the time average between tⁿ and tⁿ⁺¹. ${\tilde{\mathbf{q}}}_{i,j,k}^n$ denotes the vector of conservative variables after flux integration.

Except for the pressure gradient term of the momentum equation, which was treated during the calculation of numerical fluxes in ATHENA++, all source terms in the continuity, momentum, and energy equations can be added explicitly with an operator-splitting method, following the standard ATHENA++ approach. However, examining the source term vector in Eq. (15), all of the components in the source term, S, can be integrated analytically. Thus, instead of following the standard explicit approach, we incorporated the source terms by updating the conservative variables with the analytical solutions at tⁿ⁺¹; namely, $${\mathbf{q}}_{i,j,k}^{n+1}=\psi \left({\tilde{\mathbf{q}}}_{i,j,k}^n\right),$$(17)

where $$\psi (\mathbf{q})=\left(\begin{array}{c}{\psi }_1\\ {\psi }_2\\ {\psi }_3\\ {\psi }_4\\ {\psi }_5\end{array}\right)=\left(\begin{array}{c}\frac{J^n}{J^{n+1}}{q}_1\\ {\left(\frac{{\mathcal{R}}^n}{{\mathcal{R}}^{n+1}}\right)}^2\frac{J^n}{J^{n+1}}{q}_2\\ {\left(\frac{{\mathcal{R}}^n}{{\mathcal{R}}^{n+1}}\right)}^2\frac{J^n}{J^{n+1}}{q}_3\\ {\left(\frac{L_z^n}{L_2^{n+1}}\right)}^2\frac{J^n}{J^{n+1}}{q}_4\\ {\psi }_1{\mathcal{E}}_{k,\text{ rell }}^{n+1}+q_1\epsilon {\left(\frac{J^n}{J^{n+1}}\right)}^{\gamma }\end{array}\right),$$(18)

with $J^n\equiv {\mathcal{R}}^n{\mathcal{R}}^n{L}_z^n$ the Jacobian determinant in the local coordinate system, and $${\mathcal{E}}_{k,\text{ rell }}^{n+1}=\frac{1}{2{\psi }_1^2}{\left(L_z^{n+1}\right)}^2{\psi }_4^2+\frac{1}{2{\psi }_1^2}{\left({\mathcal{R}}^{n+1}\right)}^2({\psi }_2^2+{\psi }_3^2)$$(19)

the kinetic energy density of relative motion at tⁿ⁺¹. This approach of analytical integration of the source term significantly enhances the numerical accuracy, and ensures that the conservation laws are fully obeyed (see Sec. 4.3). The method we used to integrate the source term is first-order in time. While higher-order methods have been explored, they are not as efficient within the current framework. We further discuss the impact of the treatment of the source term on numerical accuracy in Sec. 4.4.

The flux integration (or the time-averaged numerical fluxes) mentioned above were computed depending on different time-integrating methods (see Sect. 3.2.3 of Stone et al. 2020, for details). To explain the core of the algorithm, we present hereafter the flux calculation in the x direction only as an example, dropping the subscripts j and k for simplicity. For a time integrator with a stage number of N_stage in each time step, for each stage $s=0,\dots N_{\text{stage}-1},{\mathbf{q}}_i^n$ was updated with fluxes ${\mathbf{f}}_{i-1/2}^s$ and ${\mathbf{f}}_{i-1/2}^s$, which were computed with a Riemann solver. In ATHENA++, various solvers were implemented, including Roe’s linearized solver (Roe 1981), as well as the Harten–Lax–van Leer (HLL) approximated Riemann solvers (Harten et al. 1983). These approximate Riemann solvers provide an estimate of the fluxes. As was discussed above, pressure gradient terms in the momentum equation were integrated during the calculation of fluxes in ATHENA++. Since the corresponding terms in the local model (the right-hand side of Eqs. (6)–(8)) can be interpreted as an anisotropic sound speed, we modified the Riemann solver, which incorporates the anisotropic effective sound speed by rescaling the wave speed during flux calculation. We applied this modification in Roe’s solver, the Harten–Lax–van Leer–Einfeldt (HLLE) solver, and the Harten–Lax–van Leer–Contact (HLLC) solver, the latter being recommended for adiabatic hydrodynamics. The basic introduction of these Riemann solvers can be found in Toro (1999), so will not be introduced in detail.

For the flux, ${\mathbf{f}}_{i-1/2}^s$, at the cell interface, x_i−1/2, the Roe’s solver and HLLE solver provide the estimated flux, ${\mathcal{F}}_{i-1/2}^{\text{Roe }}$ and ${\mathcal{F}}_{i-1/2}^{\text{HLLE }}$, separately. The Roe’s flux is $${\mathcal{F}}_{i-1/2}^{\mathrm{R}\mathrm{o}\mathrm{e}}=\frac{1}{2}({\mathbf{f}}_{L,i-1/2}+{\mathbf{f}}_{R,i-1/2}+{\displaystyle \sum _{\alpha }}{a}^{\alpha }\left|{\lambda }^{\alpha }\right|{\mathcal{R}}^{\alpha }).$$(20)

The HLLE flux at the interface, x_i−1/2, is $${\mathcal{F}}_{i-1/2}^{\mathrm{H}\mathrm{L}\mathrm{L}\mathrm{E}}=\frac{b^{+}{\mathbf{f}}_{L,i-1/2}-b^{-}{\mathbf{f}}_{R,i-1/2}}{b^{+}-b^{-}}+\frac{b^{+}{b}^{-}}{b^{+}-b^{-}}({\mathbf{q}}_i-{\mathbf{q}}_{i-1}).$$(21)

In both cases, the subscripts L and R denote the left and right states of each variable at the interface, and f_L,i−1/2 = f(q_L,i−1/2) and f_R,i−1/2 = f (q_R,i−1/2) are the fluxes evaluated using the left and right states of the conserved variables.

For Roe’s solver, λ^α denotes the α-th eigenvalue of the Roe’s matrix in the conserved variables (see e.g. Sect. 4.3.2 and Appendix B of Stone et al. 2008), and $$a^{\alpha }={\mathcal{L}}^{\alpha }\cdot ({\mathbf{q}}_{L,i-1/2}-{\mathbf{q}}_{R,i-1/2}),$$(22)

where the ℒ^α and ℛ^α are the αth left and right eigenvectors of Roe’s matrix, which corresponds to λ^α. ℒ^α and ℛ^α are also rows and columns of the left and right eigenmatrices, ℒ and ℛ. For isothermal hydrodynamics, the eigenvalues of Roe’s matrix are $$\lambda =(v_x-C,v_x,v_x,v_x+C),$$(23)

where C is the (effective) sound speed. The left and right eigenmatrices are $$\mathcal{L}=\left[\begin{array}{cccc}(1+v_x/C)/2& -1/(2C)& 0& 0\\ -v_y& 0& 1& 0\\ -v_z& 0& 0& 1\\ (1-v_x/C)/2& 1/(2C)& 0& 0\end{array}\right]$$(24)

and $$\mathcal{R}=\left[\begin{array}{cccc}1& 0& 0& 1\\ v_x-C& 0& 0& v_x+C\\ v_y& 1& 0& v_y\\ v_z& 0& 1& v_z\end{array}\right].$$(25)

Here, we rescaled the effective sound speed with a scaling factor, l, so that C = c_s/l, with c_s the isothermal or adiabatic sound speed. In the x and y directions, we chose l = ℛ so that C = c_s,H; in the z direction, l = L_z so we have C = c_s,V. We note that when l = 1, we returned to the original eigenvalues and eigenvectors given in Stone et al. (2008). For adiabatic hydrodynamics, the eigenvalues and the eigenmatrice are a bit more complicated than the isothermal case as the (modified) energy equation is also solved. The corresponding quantities are shown in Appendix A.

Similarly, for the HLLE solver, $$\begin{array}{rl}{b}^{+}& =max[max({\lambda }^M,v_{x,\mathrm{R}}+C_{\mathrm{R}}),0],\\ b^{-}& =min[min({\lambda }^0,v_{x,\mathrm{L}}-C_{\mathrm{L}}),0].\end{array}$$(26)

Here, λ^M and λ⁰ denote the maximum and minimum eigenvalues of Roe’s matrix given above, υ_x,L and υ_x,R are the velocity component normal to the interface in the left and right states, respectvely, and C_L = c_s,L/l, C_R = c_s,r./l are the effective sound speed computed from the left and right states. For the HLLC solver, the middle contact wave is also included. To compute the fastest wave speed including the middle wave, the rescaled sound speeds, C_L and C_R, were also used to estimate the pressure in the intermediate states. In the middle wave speed estimation processes, the pressure terms, p*, that originate from the momentum equations were replaced with p*/l².

The spatial reconstruction makes use of the left and right eigenmatrices of the Roe’s matrix in the primitive variables, and was modified accordingly. For isothermal hydrodynamics, the eigenmatrices are identical to the original ones presented in Appendix B of Stone et al. (2008), with the isothermal sound speed replaced by the effective sound speed, C, in our case, as was illustrated previously. For adiabatic hydrodynamics, the eigenmatrices are again in a different form, and are shown in Appendix A. We implemented this modification for both piecewise linear (PLM) and piecewise parabolic (PPM) spatial reconstructions.

Rescaling the sound speed implies the modification of the Courant–Friedrichs–Lewy (CFL) condition. We further considered the effect of the time-varying sound speed on the stability condition, by applying the leapfrog finite difference equation (FDE) to a harmonic oscillator with a time-varying frequency. Our new CFL condition is thus given by $$\delta t=C_{\circ }\cdot min[\frac{\delta xv\left({\beta }_x\right)}{\left|v_x\right|+c_{\mathrm{s},\mathrm{H}}},\frac{\delta yv\left({\beta }_y\right)}{\left|v_y\right|+c_{\mathrm{s},\mathrm{H}}},\frac{\delta zv\left({\beta }_z\right)}{\left|v_z\right|+c_{\mathrm{s},\mathrm{V}}}]\text{, }$$(27)

where C_∘ is the CFL number², and the original physical sound speed has been replaced by the effective sound speeds, c_s,H and c_s,V. The effect of the time-varying sound speed is reflected in υ(β_x), υ(β_y), and υ(β_z), with β_x = β_y = − U₀/c_s, and β_z = − U_r0 L_z/c_s. Both β and υ(β) are dimensionless. For a collapsing system (β > 0), which is the focus of this work, we have $$v(\beta )=\frac{\sqrt{1+8\beta }-1}{4\beta }.$$(28)

For a static system, β = 0, we have υ(β → 0) = 1, which goes back to the original CFL condition, but with rescaled sound speeds. The case of the expanding system (β < 0), as well as the derivation of the new CFL condition, are described in Appendix B.

This local model is physically meaningful with periodic boundary conditions. Further extension of the shear-periodic boundary condition (local model with weak rotation) will be considered in future work.

4 Numerical tests

In this section, we present benchmark tests of the local box for nonlinear solutions (Sec. 4.1) and linear waves (Sec. 4.2). Fiducially, we used Roe’s solver with the van Leer integrator (VL2) and PPM spatial reconstruction for our showcase tests, with an isothermal EoS. Other parameters of simulation setups, including simulation box sizes and resolutions, are described in each subsection (Sects. 4.1, 4.2, and 4.3) for each individual tests. Additionally, we present convergence studies in Sec. 4.4, examining both Roe and HLL solvers, considering various time integrators including VL2, second- and third-order Runge–Kutta (RK2, RK3) integrators, with both PPM and PLM reconstructions, and both isothermal and adiabatic EoSs (with γ = 1.4) tested. The CFL number was set to C_∘ = 0.4 for shear flow tests (Sec. 4.1) and C_∘ = 0.2 for wave tests (Sec. 4.2). We set the code unit with an initial background (unperturbed) gas density of ρ₀ = 1 and a physical sound speed of c_s0 = 1. The length unit of the code, L₀ = 1, was set by the unit of the local coordinate system defined in Eq. (1) (see also the right panel of Fig. 1). We also defined the initial values of the length scales in the local model, R₀ ≡ ℛ(0), L_z0 ≡ L_z(0), for convenience. For each numerical test, we compared our numerical results with analytical solutions to validate the accuracy and robustness of our implementation. We briefly summarize the analytical solutions corresponding to the tests in each subsections separately. Details of the analytical solutions are presented in LL23.

4.1 Nonlinear solutions

4.1.1 Horizontal shear flows

Consider a local model with a spatially homogeneous density, ρ = ρ(t), and a purely horizontal velocity field, $$v_x=A(t)v_{x0}(1+\cos \left(k_{z}z\right)),$$(29) $$v_y=A(t)v_{y0}(1+\sin \left(k_{z}z\right)),$$(30) $$v_z=0\text{. }$$(31)

This is a horizontal shear flow that varies only in the vertical direction. Without loss of generality, we can set A(0) = 1, and the solution to this local model is given with the evolution of the flow amplitude, A, $$A={(\mathcal{R}(t)/R_0)}^{-2},$$(32)

yielding the following velocity field at time, t, $$v_x=v_{x0}(1+\cos \left(k_{z}z\right)){(\mathcal{R}(t)/R_0)}^{-2},$$(33) $$v_y=v_{y0}(1+\sin \left(k_{z}z\right)){(\mathcal{R}(t)/R_0)}^{-2}.$$(34)

In global geometry, this growth in the horizontal flow velocity is a byproduct of the conservation of angular momentum during the spherical collapse.

We tested this solution by setting a local model of uniform collapse, with R₀ = 10.0L₀, U₀ = −1.0c_s0, L_z0 = 10.0L₀, and U_R0 = 0.0. The initial flow structure was set with υ_x0 = 0.025c_s0, υ_y0 = 0.0125c_s0, and k_z = 4π. The simulation box size was set with L_x,sim = L_y,sim = 0.5L₀, L_z,sim = 1.0L₀, with a resolution of 128 cells per unit length. We note that our choice of simulation box size quickly becomes too large to be physically meaningful; however, this is irrelevant to numerically testing the local model. The initial distribution of this velocity field is shown in the upper panel of Fig. 2, presenting a 3D visualization of the simulation box. Arrows mark the direction of the velocity vectors, while colors mark the velocity magnitude.

The middle panel of Fig. 2 present the simulation result with the peak shear velocity as a function of ℛ, where the peak shear velocity is measured by the maximum value of the horizontally averaged flow velocity. As was expected, the shear flow is amplified as the collapse proceeds, in excellent agreement with the analytical solution. We further defined the relative error with L2 norm of each variable, $\mathcal{X}$, $$\mathrm{err}(\mathcal{X})=\frac{1}{N}\sqrt{{\displaystyle \sum _i}{\left(\frac{{\mathcal{X}}_i-{\mathcal{X}}_{i,\text{ exact }}}{{\mathcal{X}}_{\text{max,exact }}}\right)}^2},$$(35)

where N is the total cell number, i denotes each cell, and $\mathcal{X}$_i and $\mathcal{X}$_i,exact are the numerical result and analytical solution for $\mathcal{X}$, respectvely. We used the maximum absolute value of the analytical solution over the domain, $\mathcal{X}$_max,exact = max(|$\mathcal{X}$_i,exact|), to normalize the relative error, in order to avoid meaninglessly large relative errors when the analytical solution is zero or a small number. The time evolution of err(υ_x) is displayed in the lower panel of Fig. 2, showing that our result matches the analytical solution to machine precision.

4.1.2 Elevator flows

We considered a second class of nonlinear flows in the local model, with a spatially homogeneous density, ρ = ρ(t), and a purely vertical velocity field, $$v_x=v_y=0.$$(36) $$v_z=B(t)v_{z0}(1+\cos \left(k_{x}x\right)).$$(37)

This vertical “elevator flow” has a solution with a temporally evolving flow amplitude, $$B={\left(L_z/L_{z0}\right)}^{-2},$$(38)

if without loss of generality we set B(0) = 1. This leads to a flow velocity of $$v_z=v_{z0}(1+\cos \left(k_{x}x\right)){\left(L_z/R_0\right)}^{-2}.$$(39)

These elevator flows can only occur with vertically periodic boundary condition, and could be an artifact of periodic boundaries.

To test the solution of these elevator flows, we set the local model with R₀ = 10.0L₀, U₀ = 0.0, L_z0 = 10.0L₀, and U_R0 = −0.1c_s0L₀. The flow structure was initialized with υ_z0 = 0.05c_s0 and k_x = 8π. We chose a 2D simulation box with L_x,sim = 0.5L₀ and L_z,sim = 1.0L₀, and a resolution of 256 cells per unit length. This initial flow structure is displayed in the upper panel of Fig. 3.

The middle panel of Fig. 3 presents the simulation result of the peak flow velocity as a function of L_z. The peak flow velocity was obtained by the maximum value of the vertically averaged flow velocity. The elevator flow grows with time, closely matching the analytical expectation. We further measured the relative error, err(υ_z), defined in Eq. (35). The lower panel of Fig. 3 of the time evolution of err(υ_z) shows that our result matches the analytical solution to machine precision.

Fig. 2 Test case of horizontal shear flows. Upper panel: 3D visualization of the initial spatial distribution of the shear flow in our simulation. Arrows mark the direction of shear flow, and the total flow velocity magnitude is mapped to colors (see the color bar). Middle panels: solid red line shows υx/υx0 vs. R, i.e., normalized shear velocity as a function of global radius of the local box, showing growth of shear velocity as collapse proceeds, closely matches the analytical expectation (dashed black lines). The result of υy is identical to that of υx and is not shown here. Lower panel: time evolution of relative error err(υx).

4.1.3 Diagonal flows

We next considered a class of local models with a spatially homogeneous density of ρ = ρ(t) and nonzero horizontal and vertical velocities. The solution for these diagonal flows combines A = (ℛ/R₀)⁻² of the horizontal shear flows and B = (L_z/R₀)⁻² of the elevator flows such that $$v_x={\left(\mathcal{R}/R_0\right)}^{-2}{v}_{x0}(1+\cos \left(k_{y}y\right)),$$(40) $$v_z={\left(L_z/L_{z0}\right)}^{-2}{v}_{z0}(1+\cos \left(k_{y}y\right)),$$(41)

where, without loss of generality, we have aligned the diagonal flow with the x axis. In the global picture, these diagonal flows can be interpreted as a spiral with the pitch angle evolving over time as the collapse proceeds, since the horizontal and vertical components of the flow velocity change independently.

We tested this solution by setting up a local model with R0 = 10.0L₀, U₀ = −1.0c_s0, L_z0 = 10.0L₀, and U_R0 = −0.1c_s0/L₀, with an initial flow structure of υ_x0 = 0.025c_s0, υ_z0 = 0.05c_s0, and k_y = 8π. The simulation box was set with a box size of L_x,sim = L_y,sim = 0.5L₀, L_z,sim = 1.0L₀, and a resolution of 128 cells per unit length. Figure 4 presents the flow structure at the beginning of the simulation (ℛ = R₀), as well as later on at ℛ = 0.3R₀. Variation in the flow direction is observed as the collapse proceeds. As U₀ > U_R0L_z, the horizontal flow component grows faster than the vertical component, causing the flow orientation to evolve toward the horizontal direction. The time evolution of flow amplitude closely matches the analytical solution in both the horizontal and vertical directions. The results are largely similar to those shown in Sects. 4.1.1 and 4.1.2, and are not presented here.

4.2 Linear perturbation (waves)

We examined the properties of linear perturbations in the local model (for details of the linear wave solutions refer to Sect. 4 and Appendix B of LL23; these are briefly summarized here). The density and velocity perturbations, δρ and δv, of a diagonally propagating sound wave can be written as $$\delta v_x=-\frac{2c_s}{{\mathcal{R}}^2}{k}_x\mathfrak{R}[X\exp (i\mathbf{k}\cdot \mathbf{x})],$$(42) $$\delta v_z=-\frac{2c_s}{L_z^2}{k}_z\mathrm{\Re }[X\exp (i\mathbf{k}\cdot \mathbf{x})],$$(43) $$\delta \rho =\frac{2\rho }{c_{\mathrm{s}}}\mathfrak{J}[\mathrm{\Pi }\text{\hspace{0.17em}exp}(i\mathbf{k}\cdot \mathbf{x})],$$(44)

where the (complex) wave amplitudes (X, Π) fulfill $$\dot{X}=\mathrm{\Pi },$$(45) $$\dot{\mathrm{\Pi }}=-{\omega }^{2}X.$$(46)

Here, we have introduced the sound wave frequency $$\omega =c_{\mathrm{s}}|\mathbf{k}|=c_{\mathrm{s}}{({\mathcal{R}}^{-2}{k}_x^2+L_z^2{k}_z^2)}^{1/2}.$$(47)

Without loss of generality, we aligned the wave propagating direction with the x − z plane so that δυ_y = 0.

The behavior of the wave is characterised by two regimes, depending on the value of the ratio between the wave period and the evolution timescale of the background flow, $t_{\mathrm{b}\mathrm{g}}\sim min(\mathcal{R}/\dot{\mathcal{R}},L_z/{\dot{L}}_z)$. For a slowly varying background flow with ωt_bg ≫ 1, the wave is described well by a Wentzel–Kramers–Brillouin (WKB) approximation. The WKB approximation for the density and velocity perturbations, δρ and δν, of a diagonally propagating sound wave is $$\delta v_x\approx -\frac{2c_{\mathrm{s}}}{{\mathcal{R}}^2}{\omega }^{-1/2}{k}_x\mathfrak{R}[X_{\pm }\exp (i\mathbf{k}\cdot \mathbf{x}\pm i\int \omega \mathrm{d}t)],$$(48) $$\delta v_z\approx -\frac{2c_{\mathrm{s}}}{L_z^2}{\omega }^{-1/2}{k}_z\mathrm{\Re }[X_{\pm }\exp (i\mathbf{k}\cdot \mathbf{x}\pm i\int \omega \mathrm{d}t)],$$(49) $$\delta \rho \approx \frac{2\rho }{c_{\mathrm{s}}}{\omega }^{1/2}\mathfrak{R}[X_{\pm }\exp (i\mathbf{k}\cdot \mathbf{x}\pm i\int \omega \mathrm{d}t)],$$(50)

where X_± is some (complex) constant.

For a more general test, the exact solution for the sound waves can be obtained for a given collapse profile. One defines a length scale of reference L = L(t), and adopts $$L=L_0\frac{\mathcal{R}}{{\mathcal{R}}_0}\left(\frac{k_x^2+b_0^{-2}{k}_z^2}{k_x^2+b^{-2}{k}_z^2}\right),$$(51)

where b = L_z/ℛ is the aspect ratio of the local reference frame, and the subscription “0” denotes the initial values. This length scale of reference can be equal to the global radius of the cloud, L = ℛ, which can be achieved by adopting L₀ = ℛ₀, and which can have either b = b₀ (constant aspect ratio) or k_z = 0 (purely horizontal waves). If we consider a collapse profile in the form of L = L₀ (1 − t/t_c)^β, where t_c denotes the collapse timescale, this profile includes two typical types of collapse profiles; namely, uniform collapse (L = L₀ + U₀t) where β = 1, and a free-fall profile where β = 2/3. For uniform collapse, the exact solution for density and velocity is $$\delta v_x=\frac{2c_{\mathrm{s}}}{{\mathcal{R}}^2}{k}_x\sqrt{L}\mathrm{\Re }[A_{\pm }{L}^{\pm \mathrm{\Lambda }}\exp (i\mathbf{k}\cdot \mathbf{x})],$$(52) $$\delta v_z=\frac{2c_{\mathrm{s}}}{L_z^2}{k}_z\sqrt{L}\mathrm{\Re }[A_{\pm }{L}^{\pm \mathrm{\Lambda }}\exp (i\mathbf{k}\cdot \mathbf{x})],$$(53) $$\delta \rho =\frac{2\rho L_0}{c_{\mathrm{s}}^2{t}_c}\sqrt{L}\{\pm \mathfrak{J}[(\mathrm{\Lambda }\pm 1/2)A_{\pm }{L}^{\pm \mathrm{\Lambda }}\exp (i\mathbf{k}\cdot \mathbf{x})]\}.$$(54)

Here, we have introduced $\mathrm{\Lambda }=\frac{1}{2}\sqrt{1-4{\left(c_{\mathrm{s}}\left|{\mathbf{k}}_0\right|t_c\right)}^2}$, and |k₀| is the initial amplitude of k.

For a general power law profile with β ≠ 1, the exact solution can be written as $$\begin{array}{rl}\delta v_x=& \frac{2k_x}{{\mathcal{R}}^2}{s}^{1/[2(1-\beta )]}\{\mathfrak{R}[A\exp (i\mathbf{k}\cdot \mathbf{x})]J_{\nu }(s)\\ & +\mathfrak{R}[B\exp (i\mathbf{k}\cdot \mathbf{x})]Y_{\nu }(s)\},\end{array}$$(55) $$\begin{array}{rl}\delta v_z=& \frac{2k_z}{L_z^2}{s}^{1/[2(1-\beta )]}\{\mathrm{\Re }[A\exp (i\mathbf{k}\cdot \mathbf{x})]J_{\nu }(s)\\ & +\mathrm{\Re }[B\exp (i\mathbf{k}\cdot \mathbf{x})]Y_{\nu }(s)\},\end{array}$$(56) $$\begin{array}{rl}\delta \rho =& -\frac{2\rho }{c_{\mathrm{s}}^2}\frac{1-\beta }{t_0}{s}^{\frac{1-2\beta }{2(1-\beta )}}\mathrm{sgn}(s)\{\mathfrak{I}[A\exp (i\mathbf{k}\cdot \mathbf{x})]J_{\nu -1}(s)\\ & +\mathfrak{J}[B\exp (i\mathbf{k}\cdot \mathbf{x})]Y_{\nu -1}(s)\},\end{array}$$(57)

where we have defined $$s={\left(\frac{t_c-t}{t_0}\right)}^{1-\beta },t_0={\left[\frac{|1-\beta |}{c_{\mathrm{s}}\left|{\mathbf{k}}_0\right|t_c^{\beta }}\right]}^{1/(1-\beta )},$$(58)

and J_ν and Y_ν correspond to Bessel functions of the first and second kinds, respectvely, with an order of ν = 1/ (2|1 − β|).

Additionally, given the relation $t_{\mathrm{b}\mathrm{g}}=L/|\dot{L}|(\beta \ne 1)$, the collapsing system can evolve from the WKB regime (ωt_bg ≫ 1) to the freeze-out regime (ωt_bg ≲ 1), where the wave has a longer period than the background flow timescale, and can thus show a roughly static pattern. The corresponding length scale at freeze-out can be obtained by setting ωt_bg = 1: $$L_{\text{freezeout }}=L_0{\left(\frac{t_0}{t_c}\right)}^{\beta }{\left(\frac{|\beta |}{{\omega }_0{t}_0}\right)}^{\beta /(1-\beta )}.$$(59)

For a diagonally propagating sound wave, if the aspect ratio b ≠ 1, the sound wave generates a shear flow perpendicular to the wave propagation direction in the local reference frame. This shear flow can be characterized by $$k_z\delta v_x-k_x\delta v_z=-\frac{2c_{\mathrm{s}}}{{\mathcal{R}}^2}{k}_x{k}_z(1-b^{-2})\mathfrak{R}[X\exp (i\mathbf{k}\cdot \mathbf{x})].$$(60)

For a constant aspect ratio, b, this shear flow is not a physical phenomenon, but rather simply a consequence of the choice of a spatially anisotropic reference frame. However, with a varying aspect ratio during the collapse, the shear flow generated by the sound wave has a physical origin, as it cannot be removed by a time-independent rescaling of the anisotropic reference frame.

These aforementioned wave properties in the local model deserve careful testing numerically. Here, we tested the solutions for a horizontal sound wave with both a uniform and a general power law collapse profile separately, as well as the solution for diagonal sound waves.

Fig. 3 Same as Fig. 2, but for the test case of elevator flows. The vertical flow velocity grows as the simulation box collapses, closely matching the analytical expectation.

Fig. 4 3D visualization of the spatial distribution of the diagonal flow test. Arrows mark the direction of the flow, and the total flow velocity magnitude is mapped to colors (see the color bar). Left panels are the initial setup (ℛ = ℛ0); right panels are at ℛ = 0.3R0. Flow direction varies as the collapse proceeds.

4.2.1 Horizontal sound waves with uniform collapse

We first considered a purely horizontally propagated traveling wave with a uniform collapse profile. This can also provide a comprehensive view of the wave behavior in the local box. We set up the following perturbation: $$\delta v_{x0}=A_0{R}_0^{-2}{k}_x\sin \left(k_{x}x\right),$$(61) $$\delta v_{y0}=\delta v_{z0}=0,$$(62) $$\delta {\rho }_0=A_0\frac{{\rho }_0}{c_{\mathrm{s},\mathrm{H}0}}{R}_0^{-2}{k}_x\sin \left(k_{x}x\right),$$(63)

where A₀ is some (real) constant, and c_s,H0 = c_s0/R₀ is the initial value of the effective horizontal sound speed defined in Sec. 2. This perturbation corresponds to wave amplitudes of (X, Π) = $(\frac{i{A}_0}{2c_{\mathrm{s}}},\frac{A_0{k}_x}{2R_0})$.

For the numerical test of the local model, we chose a collapse profile with R₀ = 10.0L₀, U₀ = 1.0c_s0, L_z0 = 1.0L₀, and U_R0 = 0, with a perturbation amplitude of A₀ = 10⁻⁶, k_x = 16π. The simulation box size was set as L_x,sim = 0.5L₀, L_z,sim = 1.0L₀, with a resolution of 256 cells per unit length.

The relative amplitude of gas density perturbation, δρ/ρ, was measured with respect to the background gas density, ρ, at each time point, which was obtained by averaging the gas density over the domain. The space–time plot of the vertically averaged value of δρ/ρ presented in the upper panel of Fig. 5 shows the evolution of the wave pattern. The wave propagates rightward with a time-varying effective sound speed, c_s,H = c_s0/ℛ. As the collapse proceeds, the wave propagation is accelerated as ℛ decreases, matching the overlaid curves of theoretical phase shift profiles over time. This time-varying effective sound speed is not physical in the global picture, but is instead due to our choice of a time-dependent reference frame in the local model.

Our choice of collapse profile results in the wave operating in the WKB regime. It is thus appropriate to compare the evolution of the wave amplitude to that expected from the WKB approximation solution. The lower panel of Fig. 5 shows this comparison of the relative amplitude of density perturbation, A_ρ; = max (δρ/ρ). The numerical result of the A_ρ evolution agrees well with the WKB approximation solution, Aρ = A_ρ,0ℛ^−1/2.

Similarly, fora purely vertical sound wave, the corresponding length scale becomes L_z instead of ℛ. We have also numerically tested the vertical sound wave and found that the effective vertical sound speed evolves following c_s,H = c_s0/L_z, and the relative amplitude of δρ follows $A_{\rho }=A_{\rho ,0}{L}_z^{-1/2}$, which aligns with the theoretical expectations.

4.2.2 Horizontal waves with nonuniform collapse

We next tested the horizontal sound waves with a more general power law collapse profile. In order to compare with the exact solution Eqs. (55)–(57), we considered a collapse profile with β = 1/2 and t_c = 1, i.e., ℛ = R₀(1 − t)^1/2, and correspondingly initiated the perturbation of velocity and density as $$\delta v_{x0}=A_0\frac{k_x}{R_0{}^2}s{Y}_1\left(s_0\right)\cos \left(k_{x}x\right),$$(64) $$\delta v_{y0}=\delta v_{z0}=0,$$(65) $$\delta {\rho }_0=A_0\frac{{\rho }_0}{2t_0{c}_{\mathrm{s}0}^2}{Y}_0\left(s_0\right)\sin \left(k_{x}x\right),$$(66)

where according to Eq. (58), s₀ = s(0) = (t_c/t₀)^1−β, and t₀ can again be obtained from Eq. (58). A₀ is some real constant, and this perturbation corresponds to A = 0, B = A₀/2 in Eqs. (55)–(57). For the numerical test, we set R₀ = 1.0L₀, L_z0 = 1.0L₀, U_R0 = 0, with a perturbation of A₀ = 10⁻⁵, k_x = 2π, and a simulation box size of L_x,sim = 2.0L₀, L_z,sim = 1.0L₀, with a resolution of 256 cells per unit length.

Similar to Sect. 4.2.1, the relative amplitude of the gas density perturbation, δρ/ρ, was measured at each time point, where ρ was obtained by averaging the gas density over the domain. The space–time plot of the vertically averaged relative gas density, δρ/ρ, displayed in the upper panel of Fig. 6 shows that this perturbation leads to a horizontal standing wave. This wave exhibits phase reversal with each half wave period. We further measured the relative amplitude of the density perturbation by performing a Fourier transform to the wave pattern, and identified the amplitude of the dominant frequency component. The time evolution of the relative amplitude, δρ/ρ, is shown in the lower panel of Fig. 6, with the length scale corresponding to the freeze-out regime, L_freezeout (see Eq. (59)), marked with a vertical dotted line. The result is in good agreement with the analytical solution in both the WKB regime (ℛ > L_freezeout) and the freeze-out regime (ℛ < L_freezeout). The L2 norm of the relative error for the relative amplitude, err(δρ/ρ), defined in Eq. (35) is on the order of 10⁻⁵ when the absolute value of the amplitude solution is much larger than 0. When the solution of the amplitude is close to 0 (i.e., close to the period of phase reversal), the relative error, err(δρ/ρ), can be up to the order of 10⁻², partly due to the division of a small analytical solution value, and partly due to inherent limitations and numerical inaccuracies that can arise near the phase reversal period. The slightly diffused wave pattern during the phase reversal can also be seen in Fig. 6.

Fig. 5 Test case of purely horizontal sound wave with uniform collapse. Upper panel: space–time plot for vertically averaged relative density variation. Dashed lines show the theoretical expectation for phase evolution, in good agreement with the numerical results. Lower panel: relative amplitude of perturbation, δρ/ρ, as a function of R from the simulation (solid red line). The wave operates under the WKB regime, and is closely consistent with the theoretical expectation from the WKB approximation solution (dashed black line). The wiggles in the simulation results are not present in the WKB solution, but are consistent with the exact solution presented in Eq. (54).

Fig. 6 Test case of purely horizontal sound wave with a nonuniform collapse profile. Upper panel: space-time plot for vertically averaged relative density variation. Lower panel: relative amplitude of perturbation, δρ/ρ, as a function of R from the simulation (solid red line), consistent with the analytical solution (dashed black line). The vertical dotted line marks the length scale, Lfreezeout, corresponding to the freeze-out regime.

4.2.3 Diagonal sound waves

Diagonally propagating sound waves may introduce numerical complexities, as they generate shear flows perpendicular to the direction of wave propagation. This is illustrated in the upper panel of Fig. 4, with the wave pattern characterized by the density.

We tested diagonal sound waves by generating density and velocity perturbations of $$\delta v_{x0}=A_0{R}_0^{-2}{k}_x\sin (k_{x}x+k_{z}z),$$(67) $$\delta v_{y0}=0,$$(68) $$\delta v_{z0}=A_0{L}_{z0}^{-2}{k}_z\sin (k_{x}x+k_{z}z),$$(69) $$\delta {\rho }_0=A_0\frac{{\rho }_0}{c_{\mathrm{s}0}}\sqrt{k_x^2{R}_0^{-2}+k_z^2{L}_{z0}^{-2}}\sin (k_{x}x+k_{z}z),$$(70)

where A₀ is some (real) constant. This perturbation corresponds to an initial wave amplitude of (X₀, Π₀) = $(\frac{i{A}_0}{2c_{s0}},\frac{A_0}{2}\sqrt{k_x^2{R}_0^{-2}+k_z^2{L}_{z0}^{-2}})$. Recalling Eqs. (45), (46), and (60), this pure diagonal sound wave generates a shear flow that can be characterized by $$\delta v_{\text{shear }}=k_z\delta v_x-k_x\delta v_z$$(71) $$=-\frac{2c_{\mathrm{s}}}{{\mathcal{R}}^2}{k}_x{k}_z(1-b^{-2})\mathrm{Re}[X\exp (i\mathbf{k}\cdot \mathbf{x})],$$(72)

and the dynamics of the system of (X, Π) can be obtained from Eqs. (45) and (46). Thus, the amplitude X can be obtained by solving the equation $$\ddot{X}=-{\omega }^2(t)X,$$(73)

where the wave frequency, ω, is given in Eq. (47).

We numerically tested the diagonal sound wave by setting R₀ = 1.0L₀, U₀ = 1.0c_s0, L_z0 = 1.0L₀, and U_R0 = 0, with a perturbation of A₀ = 10⁻⁵, k_x = 2π, k_z = 4π, and a simulation box of L_x,sim = 2.0L₀, L_z,sim = 1.0L₀, with a resolution of 512 cells per unit length.

In the lower panel of Fig. 7, we compare the amplitude of the shear flow generated by the diagonal sound wave in the numerical simulation with the analytical solution. The numerical result is obtained by calculating the shear flow, δυ_shear, in each grid, and performing a 2D Fourier transform of the shear flow map, followed by identifying the amplitude of the dominant frequency component. The analytical solution of the shear flow amplitude, 2c_sℛ⁻²k_xk_z (1 − b⁻²)X, was obtained by numerically solving for the wave amplitude, X, at each time step according to Eq. (73), with the explicit Runge–Kutta integration method of order 5(4) (Dormand & Prince 1980). The numerical result and analytical solution of the shear flow amplitude are in good agreement. We further measured the L2 norm of the relative error for the shear strength, δυ_shear, in each grid, and got err(δυ_shear) as defined in Eq. (35) of the order of 10⁻⁶.

Fig. 7 Test case for diagonal sound wave. Upper panel: illustration of the diagonal sound wave and the shear flow generated. Lower panel: amplitude of the shear flow generated by the wave as a function of R. The numerical result (solid red line) closely matches the analytical solution (dashed black line).

4.3 Conservation laws

The continuity equation of the local model, Eq. (5), gives the conservation law for the total mass, $$M=\iiint \rho \mathrm{d}V=\iiint \rho J\mathrm{d}x\mathrm{d}y\mathrm{d}z,$$(74)

where J ≡ ℛ² L_z is the Jacobian determinant in the local coordinate system.

The momentum equation of the local model, Eqs. (6)–(8), gives the conservation law for the total momenta, $${\mathrm{\Pi }}_i=\iiint \rho {\mathcal{V}}_i\mathrm{d}V=\iiint \rho {\mathcal{V}}_{i}J\mathrm{d}x\mathrm{d}y\mathrm{d}z,$$(75)

where ${\mathcal{V}}_i$ is the covariant velocity component in each direction, $${\mathcal{V}}_x=v_x{\mathcal{R}}^2,{\mathcal{V}}_y=v_y{\mathcal{R}}^2,{\mathcal{V}}_z=v_z{L}_z^2.$$(76)

Another conserved quantity in the local model is the kinetic helicity, H_k, of the relative flow, which is a quantity associated with the vorticity. It is defined as $$H_k=\iiint {\mathcal{V}}_i{\omega }^i\mathrm{d}V=\iiint {\mathcal{V}}_i{\omega }^{i}J\mathrm{d}x\mathrm{d}y\mathrm{d}z,$$(77)

where ωⁱ is the fluid (contravariant) vorticity, $${\omega }^i={\epsilon }^{ijk}{\mathrm{\partial }}_j{\mathcal{V}}_k,$$(78)

and ϵ^ijk is the volume element of the local model. Details and derivations of the conservation laws in the local model are discussed in Sect. 2.5 of LL23.

In order to validate these conservation laws in the local box, we performed a simulation run with a uniform background density, ρ₀ = 1, introducing white noise perturbations in the velocity field, υ_x, υ_y, and υ_z, with a relative amplitude of 10⁻⁶. We chose a collapse profile with R₀ = 10.0, U₀ = −1.0c_s0, L_z0 = 10.0, and U_R0 = −0.1, and set the simulation box size as L_x,sim = L_y,sim = 0.5, L_z,sim = 1.0, with a resolution of 32 cells per unit length. This test run is labeled as “WN.” For comparison, we also performed a simulation with the same setup for the white noise, but without collapse (i.e., U₀ = U_R0 = 0), which is labeled as “WN0”. In addition, the result of the horizontal shear flow test (see Sect. 4.1.1) was also used to examine the conservation laws.

We measured the volume-averaged density, $\overline{\rho }$, total momenta, Π, and kinetic helicity, H_k, in the simulations, and present in Fig. 8 the normalized variation in $J\overline{\rho }$, Π, and H_k as a function of time. Our results demonstrate that the modified mass conservation law ($J\overline{\rho }$) and the conservation law for total momenta (Π_i) are both satisfied to machine precision. The kinetic helicity conservation, although preserved to machine precision for the shear flow test, exhibits significant deviation for the white noise simulations. This discrepancy persists in both WN and WN0 cases, suggesting that the lack of conservation is not a consequence of the integration of the collapsing box. Instead, it is likely due to intrinsic characteristics of the ATHENA++ framework, which in certain situations is not strictly vorticity conserving, likely because the higher-order Godunovscheme generally introduces a degree of numerical diffusion (see e.g. Seligman & Laughlin 2017). In the white noise case in which vorticity is dominated by fine-scale structures, there is increased numerical dissipation at small scales because the slope limiter is expected to lead to the nonconservation of H_k. For an adiabatic EoS, the entropy conservation law is also validated; this is further presented in Appendix C.

4.4 Convergence study

We performed a convergence study with the diagonal wave (Sect. 4.2.3) test, as it shows the maximum possible error. Basic simulation setups such as wave properties and box sizes are the same as presented before. Fiducially, simulations were performed with Roe’s solver, the VL2 integrator, PPM reconstruction, and an isothermal EoS. The CFL number in the convergence study was fixed to 0.2. However, we validated the numerical stability of our code for CFL numbers up to 0.4 for VL2 and 0.8 for RK3. We further explored the effects of different Riemann solvers, time integrators, spatial reconstruction methods, and the EoS on the numerical performance of the code. We also varied the spatial resolution of our simulation in units of cells per unit length. The simulation setups for test runs that we performed for the convergence study are presented in Table 1. Figure 9 shows the time-averaged L2 relative error, err(δυ_shear), as a function of resolution.

Our results demonstrate that our code generally achieves convergence between the first and second orders for the diagonal wave test. The performances for test cases HLLE and RK3 remain generally similar to the fiducial case, while a lower order of spatial reconstruction (PLM) leads to a lower accuracy at lower resolutions. These tests with an isothermal EoS achieve second-order convergence at lower resolutions, with a tentative trend toward the first order at high resolutions. This trend is likely attributed to time-dominated errors at higher resolutions, as our source term incorporation is first-order in time. We also explored the higher-order method to integrate source terms with the VL2 integrator, which was not as efficient as our first-order method, showing a comparable or slightly higher error, and reduced stability for adiabatic flows. Therefore, we chose to use our first-order method. The source term can also be incorporated into higher-order methods with Runge–Kutta integrators. We leave the development of higher-order methods to future work.

The simulation runs with the adiabatic EoS (for both the Roe and the HLLC solvers) exhibit systematically higher (~ an order of magnitude) error levels, with first-order convergence. This pattern is likely attributed to the complexities that the adiabatic EoS introduces to the system, including the additional energy equation to be solved. Additionally, we note that our simulations may encounter stability issues in scenarios involving supersonic adiabatic flows with added perturbations (white noise). Future work will address these challenges and improve the handling of adiabatic processes within our simulations.

Fig. 8 Conservation laws for tests of horizontal shear flow (left) and of white noise in the velocity field (right). For each case, the normalized variation in $J\overline{\rho }$ (left), Πi (middle), and Hk (right) is presented as a function of time. For the horizontal shear flow test (see also Sect. 4.1.1), the modified mass conservation law ($J\overline{\rho }$) and the conservation law for total momenta (Πi) and for kinetic helicity Hk are all satisfied to machine precision (note that for all of the panels, the y axis ticks are in units of 10−13). For the case with white noise in the velocity field, the conservation laws for both mass, $J\overline{\rho }$, and total momenta, Πi, are fulfilled to machine precision as well, but not for the kinetic helicity, Hk (note that for the upper two panels, the y axis ticks are in units of 10−14). For comparison, a simulation with the same setup of white noise, but without collapse (i.e., U0 = UR0 = 0) was performed (“WN0”, dashed grey lines, see also the legend), which gives consistent results to WN. We thus anticipate that the deviation from conservation law for Hk is likely due to intrinsic characteristics of ATHENA++, which in certain situations is not strictly vorticity-conserving.

Fig. 9 Time-averaged L2 norm relative error for diagonal sound wave test, as a function of simulation resolution (cells per unit length). Overall, convergence between the first and second orders is achieved (see the dashed lines for reference).

5 Discussion

The complete series of tests presented in Sec. 4, including the challenging diagonal wave test that induces shear flows, demonstrates the robustness and stability of our ATHENA++ code implementation, confirming its capability to handle complex local model simulations effectively. ATHENA++ has also been used to incorporate similar local models; namely, the uniform expanding box for studying stellar winds (Squire et al. 2020) and the anisotropic expanding box implemented in the MHD-PIC module (Sun & Bai 2023). These previous numerical implementations have concentrated on investigating MHD characteristics and the behavior of magnetized waves. Our study focuses on a comprehensive examination of the hydrodynamic aspects of the expanding and collapsing box. Moreover, while our code is adaptable to scenarios akin to the aforementioned expanding cases, here we focus on collapsing systems. This is motivated not only by our interest in probing protostellar collapse, but also by the fact that collapsing environments – with their intrinsic amplification of density, velocity, vorticity, and energy – are potentially more abundant in instabilities. Indeed, our observations of the growth of shear flows and wave-induced shear underline the possibility of hydrodynamic instabilities during protostellar collapse. Moreover, related previous work (Robertson & Goldreich 2012) also presented an increase in turbulent velocities (“adiabatic heating”) of contracting gas. It is thus both numerically and scientifically important to conduct a thorough investigation and validation of this local collapsing box model before delving into MHD studies.

The magnetic field plays an important role in protostellar collapse, as the core must lose angular momentum during the collapse (see e.g. review of McKee & Ostriker 2007), and magnetic braking has been considered as one of the main mechanisms driving angular momentum transfer during core collapse and forming protoplanetary disks with realistic properties (e.g. Allen et al. 2003; Price & Bate 2007; Hennebelle & Fromang 2008; Commerçon et al. 2011; Masson et al. 2016; Machida et al. 2016; Wurster et al. 2016; Hennebelle et al. 2020, etc.). The extension to MHD for the theoretical model of the expanding or collapsing box will likely follow a similar approach to the MHD eccentric shearing box (Ogilvie & Barker 2014). The numerical implementation will presumably be similar to the approach described in this paper. However, magnetic fields in collapse processes may break spherical symmetry (Galli & Shu 1993; Hennebelle 2001; Galli et al. 2006), and thus may require further modifications of the local model. One should also be careful with the treatment of the background magnetic field gradient, in order to avoid spurious instabilities (Latter & Papaloizou 2017).

The local deviation from the spherical symmetry, such as that caused by the magnetic field discussed above, leads to the need for further generalization of this local box toward more realistic collapse profiles. A weak rotation in the global picture can be incorporated by simply using shear-periodic boundary conditions, whereas a strong rotation requires modification of the local model. Another factor to be considered is physical collapse profiles. As is shown in Sect. 4.2.2, a time-dependent collapse profile can transition between different regimes, including from subsonic to supersonic flows, or from WKB to freeze-out wave regimes. It is thus necessary to study the effect of various collapse profiles, such as the typical profile of Shu (1977).

Dust grains are the building blocks of planets. Recent surveys of protoplanetary disks across different stages of disk evolution (e.g. Pascucci et al. 2016; Barenfeld et al. 2016; Ansdell et al. 2018; Long et al. 2018; Cieza et al. 2019; Tychoniec et al. 2020; Tobin et al. 2020; van der Marel & Mulders 2021) show that the available solid mass to form planets is rapidly depleted during the early stage of disk evolution, likely indicating an early start of planet formation. Observational (Pagani et al. 2010; Kataoka et al. 2015; Galametz et al. 2019; Valdivia et al. 2019) and theoretical (e.g. Ormel et al. 2009; Hirashita & Yan 2009; Guillet et al. 2020; Bate 2022, etc.) works both suggest that dust grains are growing in collapsing protostellar cores, significantly prior to the protoplanetary disk phase. With further incorporation of dust, our local collapse simulation will be capable of resolving and investigating dust dynamics and their interactions with gas, especially with potential instabilities and other local phenomena, which have a rich potential to concentrate dust during the collapse. Moreover, the interaction between dust and gas drag force might give rise to additional dust-related instabilities. One example, the streaming instability (Youdin & Goodman 2005), has already been proven to be a promising avenue for planetesimal formation in protoplanetary disks, and has recently been extended to a family of instability caused by dust-gas interaction (resonance drag instability, Squire & Hopkins 2018), which could potentially operate to concentrate dust in the context of protostellar collapse. These dust concentration mechanisms in the initial collapse phase might significantly expedite the onset of planet formation, potentially starting before the complete disk formation. This aligns with the aforementioned observational indicators of early planet formation in very young stellar systems. The implementation of dust for the current model is likely similar to the method presented in this paper. However, the relative drift between the background flows of gas and dust has to be implemented in an ad hoc way, similar to the approach of the shearing box (e.g., Bai & Stone 2010).

6 Conclusion

We implemented the local collapsing box model of LL23 into ATHENA++, and performed benchmark tests for the behavior of nonlinear and linear solutions in the local model. This model features a spatially uniform time-dependent Cartesian geometry, allowing (shear-)periodic boundary conditions. The key points of this study are:

• We numerically implemented the local model by modifying the Riemann solver with a rescaled anisotropic effective sound speed and an adapted energy equation, together with operator split-source terms;

• Our benchmark test of horizontal shear flow shows that the shear is amplified as the collapse proceeds, in agreement with the analytical solution to machine precision. This corresponds to a zonal flow in the global picture, and the flow amplification reflects angular momentum conservation. Similar phenomena are verified with tests of the elevator flow and the diagonal flow, which also have corresponding physical interpretations;

• Our benchmark tests of sound waves match the theoretical predictions of wave speed and amplitude evolution over time in the local model, for both uniform and general power-law collapse profiles. Additionally, a diagonal wave generates a shear flow, which could pose numerical challenges that are well handled by our implementation. The growth and evolution of this shear flow are also rigorously tested in our simulation, with a typical relative error level of 10⁻⁴−10⁻⁶;

• The density, momentum, and vorticity in the local model increase as the collapse proceeds. The corresponding (modified) conservation laws for mass and total momenta are validated in our code;

• Our code generally achieves convergence between the first and second orders for the diagonal wave-generated shear flow test. Further improvement in performance with an adiabatic EoS is necessary for future work.

This local collapsing box implementation in ATHENA++ is not only crucial for probing the local phenomena and hydrodynamic instabilities of protostellar collapse, but also holds potential for future adaptations to include MHD processes and dust dynamics with more realistic collapse profiles. Its ability and potential to simulate local behaviors and instabilities and the interplay of various physical factors positions it as a valuable tool for advancing studies of protostellar collapse and disk formation, as well as the earliest stages of planet formation. The code will be available upon request for collaboration and will be publicly released shortly after additional benchmark tests and improvements.

Acknowledgements

We thank Francesco Lovascio, Benoît Commerçon, Xuening Bai, Xiaochen Sun, and Lile Wang for constructive discussions. We thank the anonymous referee for a thorough and helpful report. Z.X., E.L. and G.L. acknowledge the support of the European Research Council (ERC) CoG project PODCAST No. 864965.

Appendix A Adiabatic equation of state

The eigenvalues and eigenmatrices in the conserved variables with an adiabatic EoS are needed to construct the fluxes of the conserved variables in the Roe’s solver. These quantities are provided hereafter $$\lambda =(v_1-C,v_1,v_1,v_1,v_1+C),$$(A.1) $$\mathcal{L}=\left[\begin{array}{ccccc}{N}_a({\gamma }^{\mathrm{\prime }}{v}^2/2+v_{1}C)& -N_a({\gamma }^{\mathrm{\prime }}{v}_1+C)& -N_a{\gamma }^{\mathrm{\prime }}{v}_2& -N_a{\gamma }^{\mathrm{\prime }}{v}_3& N_a{\gamma }^{\mathrm{\prime }}\\ -v_2& 0& 1& 0& 0\\ -v_3& 0& 0& 1& 0\\ 1-N_a{\gamma }^{\mathrm{\prime }}{v}^2& {\gamma }^{\mathrm{\prime }}{v}_1/C^2& {\gamma }^{\mathrm{\prime }}{v}_2/C^2& {\gamma }^{\mathrm{\prime }}{v}_3/C^2& -{\gamma }^{\mathrm{\prime }}/C^2\\ N_a({\gamma }^{\mathrm{\prime }}{v}^2/2-v_{1}C)& -N_a({\gamma }^{\mathrm{\prime }}{v}_1-C)& -N_a{\gamma }^{\mathrm{\prime }}{v}_2& -N_a{\gamma }^{\mathrm{\prime }}{v}_3& N_a{\gamma }^{\mathrm{\prime }}\end{array}\right],$$(A.2) $$\mathcal{R}=\left[\begin{array}{ccccc}1& 0& 0& 1& 1\\ v_1-C& 0& 0& v_1& v_1+C\\ v_2& 1& 0& v_2& v_2\\ v_3& 0& 1& v_3& v_3\\ H-v_{1}C{L}_1^2& v_2{L}_2^2& v_3{L}_3^2& v^2/2& H+v_{1}C{L}_1^2\end{array}\right]\text{, }$$(A.3)

where N_a ≡ 1/ (2C²) and γ′ = γ − 1. In addition, in the local model we have redefined $$C=c_{\mathrm{s}}/L_1,$$(A.4) $$v^2=v_1^2{L}_1^2+v_2^2{L}_2^2+v_3^2{L}_3^2,$$(A.5)

and $$H={\mathcal{E}}_{\text{rell }}+p/\rho .$$(A.6)

The subscripts 1,2,3 in L₁,L₂,L₃ and υ₁, υ₂, υ₃ denote the corresponding quantities along the three axes of a Cartesian coordinate system, with an order that respects the right-hand orientation, i.e., in the order of (x, y, z), (y, z, x), or (z, x, y). L₁, L₂, L₃ are the characteristic length scales on the corresponding directions, which is ℛ on x and y directions, and L_z on the z direction.

The eigenvalues and eigenmatrices in the primitive variables are needed for the reconstruction algorithm. These quantities are $$\lambda =(v_1-C,v_1,v_1,v_1,v_1+C),$$(A.7) $$\mathcal{L}=\left[\begin{array}{ccccc}0& -\rho /(2C)& 0& 0& 1/\left(2C^2{L}_1^2\right)\\ 1& 0& 0& 0& -1/\left(C^2{L}_1^2\right)\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& \rho /(2C)& 0& 0& 1/\left(2C^2{L}_1^2\right)\end{array}\right],$$(A.8) $$\mathcal{R}=\left[\begin{array}{ccccc}1& 1& 0& 0& 1\\ -C/\rho & 0& 0& 0& C/\rho \\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ C^2{L}_1^2& 0& 0& 0& C^2{L}_1^2\end{array}\right].$$(A.9)

Appendix B Derivation of the CFL condition

In addition to the standard CFL condition, we further consider the effect of the time-varying sound speed on the stability condition. Consider a harmonic oscillator equation $$\ddot{y}+{\omega }^2(t)y=0,$$(B.1)

with a time-varying frequency ω (t) = ω₀ (1 + δΔt), where ω₀ is the initial frequency, and $\delta ={\dot{\omega }}_0/{\omega }_0$. Applying leapfrog FDE $$\frac{y^{n+1}-2y^n+y^{n-1}}{\mathrm{\Delta }{t}^2}+{\omega }_0^2(1+\delta \mathrm{\Delta }t{)}^2{y}^n=0,$$(B.2)

and look for solutions of the form $y^n=y_0{e}^{i\omega t}$, leading to $$e^{i\omega \mathrm{\Delta }t}-2+e^{-i\omega \mathrm{\Delta }t}+{\omega }_0^2\mathrm{\Delta }{t}^2(1+\delta \mathrm{\Delta }t{)}^2=0,$$(B.3)

which simplifies to $${\sin }^2\left(\frac{\omega \mathrm{\Delta }t}{2}\right)={\left(\frac{{\omega }_0\mathrm{\Delta }t}{2}\right)}^2(1+\delta \mathrm{\Delta }t{)}^2.$$(B.4)

This yields solutions of $$\omega =\frac{2}{\mathrm{\Delta }t}\arcsin [\pm \frac{{\omega }_0\mathrm{\Delta }t}{2}(1+\delta \mathrm{\Delta }t)],$$(B.5)

and the stability condition is thus $$|\frac{{\omega }_0\mathrm{\Delta }t}{2}(1+\delta \mathrm{\Delta }t)|\le 1.$$(B.6)

Let $x={\omega }_0\mathrm{\Delta }t>0,\beta =\frac{\delta }{{\omega }_0}$, the stability condition then becomes $$|\beta x^2+x|-2<0.$$(B.7)

For β = 0, the stability condition becomes x < 2, i.e., $\mathrm{\Delta }t<2{\omega }_0^{-1}$, which goes back to the original CFL condition. It is convenient to write the CFL condition in the form of $\mathrm{\Delta }t<2{\omega }_0^{-1}v(\beta )$, where υ(β) the additional “correction” term that reflects the time-varying sound speed, and here we have υ(β) = 1.

For β > −1/8, the stability condition eq. (B.7) yields $$x<\frac{-1+\sqrt{1+8\beta }}{2\beta },$$(B.8)

thus $$v(\beta )=\frac{-1+\sqrt{1+8\beta }}{4\beta },$$(B.9)

For β < −1/8, the stability condition eq. (B.7) yields $$x<\frac{-1-\sqrt{1-8\beta }}{2\beta },$$(B.10)

thus $$v(\beta )=\frac{-1-\sqrt{1-8\beta }}{4\beta }.$$(B.11)

Appendix C Conservation of entropy

For an adiabatic EoS in the absence of shocks, the local model leads to the conservation law for the following quantity $$\iiint \epsilon {\rho }^{2-\gamma }Jdxdydz.$$(C.1)

For a system with uniform density, this leads to conservation for the following quantity related to the entropy $$K\equiv U_{int}{M}^{-\gamma }{J}^{\gamma },$$(C.2)

where M is the total mass in the system, and $$U_{int}\equiv \iiint \rho \epsilon dxdydz,$$(C.3)

the total internal energy in the system is thus U_int/ J.

In order to validate the conservation law of entropy in the local collapsing box, we present the results for the horizontal shear flow test runs (see Sect. 4.1.1) with adiabatic EoS, using Roe and HLLC solvers, for examining the conservation law. These runs are marked with the same labels as in Table 1.

For both the Roe and HLLC solvers, K is conserved to the order of 10⁻¹⁰.

Fig. C.1 Conservation law of entropy with adiabatic EoS for tests of horizontal shear flow. For each case, the normalized variation of K are presented as a function of time. The results with the Roe’s solver are presented with red curves, while the results with the HLLC solver are presented in blue. The conservation law for entropy is satisfied to the order of 10−10.

References