Damping of MHD Turbulence in A Partially Ionized Medium

The coupling state between ions and neutrals in the interstellar medium plays a key role in the dynamics of magnetohydrodynamic (MHD) turbulence, but is challenging to study numerically. In this work, we investigate the damping of MHD turbulence in a partially ionized medium using 3D two-fluid (ions+neutrals) simulations generated with the AthenaK code. Specifically, we examine the velocity, density, and magnetic field statistics of the two-fluid MHD turbulence in different regimes of neutral-ion coupling. Our results demonstrate that when ions and neutrals are strongly coupled, the velocity statistics resemble those of single-fluid MHD turbulence. Both the velocity structures and kinetic energy spectra of ions and neutrals are similar, while their density structures can be significantly different. With an excess of small-scale sharp density fluctuations in ions, the density spectrum in ions is shallower than that of neutrals. When ions and neutrals are weakly coupled, the turbulence in ions is more severely damped due to the ion-neutral collisional friction than that in neutrals, resulting in a steep kinetic energy spectrum and density spectrum in ions compared to the Kolmogorov spectrum. We also find that the magnetic energy spectrum basically follows the shape of the kinetic energy spectrum of ions, irrespective of the coupling regime. In addition, we find large density fluctuations in ions and neutrals and thus spatially inhomogeneous ionization fractions. As a result, the neutral-ion decoupling and damping of MHD turbulence take place over a range of length scales.

Simulating two-fluid MHD turbulence is more challenging than single-fluid MHD simulations due to the high Alfvén speed of ions at low ionization fractions, which requires a much smaller time step.To address this issue, the "heavy-ion approximation" (HIA) has been adopted to accelerate explicit two-fluid MHD simulations (Oishi & Mac Low 2006;Li et al. 2006;McKee et al. 2010).This approach increases the mass of ions and reduces the ion-neutral drag coefficient  d (Draine et al. 1983;Shu 1992) accordingly.However, for simulating MHD turbulence in a weakly ionized medium, the HIA approximation may raise uncertainties (Ballester et al. 2018;Tilley & Balsara 2010).The single-fluid treatment used in, e.g., O' Sullivan & Downes (2006, 2007) for numerical modeling of MHD turbulence in a weakly ionized medium cannot fully capture the two-fluid effects in the weakly coupled regime (Tilley & Balsara 2010;Xu et al. 2016).Despite these challenges, numerical methods are crucial for testing theories of two-fluid MHD turbulence and studying ion-neutral collisional damping in an inhomogeneous medium.Three-dimensional (3D) simulations of two-fluid MHD turbulence with the RIEMANN code (Balsara 1998) have been carried out by Tilley & Balsara (2010); Meyer et al. (2014).These studies show differences in the turbulent energy spectra of ions and neutrals.The persistence of the energy cascade of Alfvén modes on scales smaller than the amplipolar diffusion scale calculated using the wave description of MHD turbulence (Burkhart et al. 2015) suggests that the damping of MHD turbulence is different from the damping of MHD waves.
In this work, we use 3D two-fluid MHD simulations to test the theoretical models developed by Xu et al. (2015Xu et al. ( , 2016) ) and study the properties of MHD turbulence in various ion-neutral coupling regimes in the presence of turbulence-induced density inhomogeneities.We perform the two-fluid MHD turbulence simulations using MHD code Athena++ (Stone et al. 2020), updated using the Kokkos framework (denoted as AthenaK).The code utilizes IMEX integrators (Pareschi & Russo 2005) to enable high-order (in time) implementation of ion-neutral drag terms, which allows for higher accuracy and stability than operator-split methods (Arzamasskiy & Stone, in prep.).To reduce the computational cost, we consider a moderately low ioniza-tion fraction.Different regimes of ion-neutral coupling are achieved by varying the numerical value of  d .To evaluate the limitations of this approach, we also carry out simulations with the same  d   but a lower ionization fraction, where   is the ion mass density.
The paper is organized as follows: in § 2, we describe the 3D numerical simulations of two-fluid MHD turbulence used in this study.In § 3, we review the recent theoretical understanding on neutral-ion decoupling and collisional damping of MHD turbulence.In § 4, we present the numerical results on the statistics of the velocity, density, and magnetic field in both ions and neutrals.The implications of the results and comparison with earlier studies are discussed in § 5, and our main findings are summarized in § 6.

Numerical setup
The 3D two-fluid simulations analyzed in this work are generated using the Athena++ implemented with Kokkos (Stone et al. 2020).
We consider the two-fluid magneto-fluid system, comprised of ions (together with electrons) and neutrals.The effects of gravity, heat conduction, ionization, and recombination are not included in the current study.The simulations solve the ideal two-fluid MHD equations, using periodic boundary conditions, IMEX3 time integration algorithm, and PPM4 spatial reconstruction method.The equations are: here  and  are the mass density and velocity of the ionized fluid (with the subscript "") and neutral fluid (with the subscript ""), respectively.We adopt an isothermal equation of state, where the sound speed   is constant.The isothermal condition applies to a medium with efficient cooling, such as molecular clouds (Tilley & Balsara 2010;Meyer et al. 2014).It only breaks down when density exceeds ∼ 10 10 cm −3 (Furuya et al. 2012).The ion-neutral collisional damping under other conditions will be studied in our future work.To drive turbulent motions in ions and neutrals, a stochastic forcing term    is utilized.Explicitly,     and     are weighted by ion and neutral densities to achieve the same injected turbulent velocities in the two fluids.
At the start of the simulation, the magnetic field and (ion and neutral) density fields are set to be uniform, with the magnetic field along the -axis.The initial ionization fraction is   =   /(  +   ), where   and   are the initial mass densities of ions and neutrals.The simulation box is divided into 480 3 cells and is uniformly staggered.

Turbulence driving
The forcing term,    , is introduced to drive the turbulence in a solenoidal manner.This is ensured by making the forcing term divergence-free.The forcing term is modeled using the stochastic Ornstein-Uhlenbeck (OU) process, which allows us to control the auto-correlation timescale,  c , of the turbulence.The auto-correlation timescale is approximately equal to  c ≈  inj /(2 A ), where  inj is the turbulence injection scale and is the Alfvén speed in the two fluids.The time step is the minimum time step allowed by the Courant-Friedrichs-Lewy stability condition for the ion and neutral fluids, respectively.
To vary the level of turbulence, we change the values of  inj .The energy injection is focused around wavenumber  = 2/ = 1 − 2 (in the unit of 2/ box , where  box is the length of simulation box) in Fourier space, where  is the length scale in real space.The turbulence is numerically dissipated at wavenumber  dis ≈ 40 − 50.We run the simulations for six eddy turnover times to ensure that the turbulence has reached a statistically stable state.
The simulation of scale-free turbulence can be characterized by the sonic Mach number,   =  inj   , and the Alfvén Mach number,   =  inj   , where  inj is the injection velocity.In this work, we fix   and   to approximate unity, ensuring that the simulations fully fall into the strong turbulence regime1 Here  dis is the turbulence dissipation scale.The critical parameters for this study are listed in Tab. 1.
Fundamental work on anisotropic incompressible MHD turbulence theory was initiated by Goldreich & Sridhar (1995) in the trans-Alfvénic regime with   ≈ 1. Goldreich & Sridhar (1995) found the "critical balance" condition, which equates the turbulence cascading time ( ⊥   ) −1 with the Alfvén wave period ( ∥   ) −1 .Here,  ⊥ is the wavevector perpendicular to the magnetic field, and   is the turbulent velocity at scale .Later studies further found that the "critical balance" condition can be valid in the strong turbulence regime when   is not unity (Lazarian & Vishniac 1999;Lazarian 2006): (i) in super-Alfvénic turbulence, where   > 1, the magnetic field's role at the injection scale,  inj , is insignificant, resulting in isotropic turbulence.However, as the turbulence cascades to smaller scales with decreasing turbulent velocity, the Alfvén speed becomes comparable to the turbulent speed at the Alfvén scale,   =  inj  −3  .This leads to the development of strong turbulence on smaller scales.(ii) in sub-Alfvénic turbulence (  < 1), the strong turbulence regime spans from the transitional scale  trans =  inj  2  to smaller scales.Turbulence within the range from  inj to  trans is termed weak Alfvénic turbulence, which is wave-like and does not obey the "critical balance".
As turbulence cascades preferentially along the direction perpendicular to the local magnetic field (Lazarian & Vishniac 1999), where the resistance to turbulent mixing of magnetic fields is the minimum, we have the scaling relations of velocity fluctuation in the strong turbulence regime: where  ⊥ is the length scale perpendicular to the local magnetic field.By using the scaling of velocity fluctuation and the "critical balance" condition, one can easily obtain the anisotropy scaling: The scale-dependent anisotropy of MHD turbulence described by the above expression indicates that the turbulent eddy has a parallel size much larger than its perpendicular size, and this anisotropy increases with the decrease of length scales.However, note that the scaledependent anisotropy can only be measured in the reference frame of local magnetic fields that percolate the turbulent eddy (Lazarian & Vishniac 1999;Cho & Vishniac 2000).

Decoupling of ions and neutrals
The interaction between ions and neutrals can be quantified by the neutral-ion collisional frequency   =  d   =  d   (  +   ) and ion-neutral collisional frequency   =  d   , respectively (Shu 1992).Neutrals start to decouple from ions when the energy cascading rate of MHD turbulence matches   .Given   ≫   in a weakly ionized medium, ions decouple from neutrals on a much smaller scale than the neutral-ion decoupling scale, so here we mainly consider neutral-ion decoupling.The coupling status between ions and neutrals in MHD turbulence can be separated into three important regimes: (i) strongly coupled regime, in which the scales are larger than the neutral-ion decoupling scale.Neutrals and ions act as singlefluid in this regime.(ii) weakly coupled regime, where the scales are The velocity maps are normalized by the mean value.The view direction is perpendicular to the mean magnetic field, which is along the vertical -direction.Bottom: Turbulent kinetic energy spectra of ions (red) and neutrals (blue).The spectra are averaged over several snapshots after turbulence reaches a statistically steady state, with the time interval equal to the largest eddy turnover time.The shadowed areas represent the variations.Magenta and green lines represent the theoretically expected neutral-ion parallel decoupling (Eq.5) and damping wavenumbers (Eq.6), respectively, for Alfvén modes of MHD turbulence.The perpendicular decoupling and damping wavenumbers are not shown because they are larger than 10 3 .smaller than the neutral-ion decoupling scale but larger than the ionneutral decoupling scale.Neutrals thus decouple from ions, but ions still couple to neutrals.(iii) Decoupled regime, in which the scales are smaller than the ion-neutral decoupling scale.Neutrals and ions in this regime are fully decoupled so if the turbulence injection happens in this regime, neutrals develop an independent hydrodynamic turbulent cascade and ions develop an MHD turbulent cascade.
In earlier linear analysis (Kulsrud & Pearce 1969), it was considered that the decoupling of neutrals from Alfvén wave oscillations at the neutral-ion decoupling wavenumber  dec, ∥ .It can be determined by equating the Alfvén wave frequency and   (Shu 1992): where the subscript "∥" means the wavevector parallel to the magnetic field.MHD turbulence was previously modeled as a collection of linear waves (Giacalone & Jokipii 1999), and  dec, ∥   =   was taken as the neutral-ion decoupling wavenumber or ambipolar diffusion wavenumber of MHD turbulence (Kulsrud & Pearce 1969;Mouschovias & Morton 1991;Hennebelle & André 2013).However, it is essential to note that MHD turbulence is a highly nonlinear phenomenon and the Alfvén wave-like motion in the strong turbulence regime with the critical balance cannot survive for more than a wave period.With the dynamically coupled turbulent mixing motion in the perpendicular direction and the wave-like motion in the parallel direction, scale-dependent anisotropy is one of its most important properties (see § 3.1).
For Alfvénic turbulence, which usually carries most of the MHD turbulence energy (Cho et al. 2002;Hu et al. 2021b), the anisotropy suggests that the neutral-ion decoupling scale is not isotropic.The parallel component of the decoupling scale can be much larger than the perpendicular component when it is significantly smaller than  inj .Taking into account the critical-balance relation between turbulent motions and wave-like motions and the anisotropy of MHD turbulence, Xu et al. (2015) derived the parallel decoupling wavenumber  dec, ∥ and perpendicular decoupling wavenumber  dec,⊥ by using the anisotropic scaling (Eq. 3) in strong MHD turbulence regimes (i.e.,  dec,⊥ >  −1  or  dec,⊥ >  −1 trans ): Here, we only consider Alfvén turbulence as it carries most of the turbulent energy (Cho & Lazarian 2003;Hu et al. 2022a).For a more in-depth discussion on neutral-ion decoupling scales of the three MHD modes (Alfvén, fast, and slow), see Xu et al. (2015Xu et al. ( , 2016)).

Neutral-ion collisional damping of MHD turbulence in a partially ionized medium
At length scales larger than  −1 dec, ∥ , ions and neutrals are perfectly coupled, and together carry the MHD turbulence.However, at length scales smaller than  −1 dec,⊥ , neutrals begin to decouple from ions, resulting in the development of their own hydrodynamic turbulent cascade, while ions continue to undergo frequent collisions with surrounding neutrals (down to the ion-neutral decoupling scale).As a result, the remaining MHD turbulence in ions is strongly affected and damped by the collisional friction exerted by neutrals, which is denoted as neutral-ion (collisional) damping.
When neutral-ion damping dominates over the damping caused by the kinematic viscosity of neutrals, the parallel damping wavenumber  dam, ∥ for Alfvénic turbulence, as derived in Xu et al. (2015) by equating the turbulent cascading rate  −1 cas =   / ⊥ and the ionneutral collisional damping rate 2  in the strong MHD turbulence regime, is given by (Xu et al. 2015(Xu et al. , 2016)): where   =   /(  +   ) is the fraction of neutrals.It holds for both sub-Alfvénic and super-Alfvénic turbulence.The perpendicular damping wavenumber  dam,⊥ can be derived from the anisotropy scaling in the strong turbulence regime: It is important to note that the  dam,⊥ is the most crucial in determining the damping of the MHD turbulent cascade because the cascade mainly happens in the direction perpendicular to the local magnetic field.In addition,  dam,⊥ is larger than  dec, ∥ (see Eq. 5), as damping of MHD turbulence takes place after neutrals decouple from ions.

Neutral-ion collisional damping of MHD turbulence
Strongly coupled regime: We present 2D slices of velocity fields for ions and neutrals in Figs. 1 and 2, taken perpendicular to the mean magnetic field at  = 240 cell.In Fig. 1, we show the cases of  d = 1 × 10 5 , 1 × 10 4 , and 1 × 10 3 , where  dec, ∥ and  dec,⊥ are expected to be larger than the numerical dissipation wavenumber We find that the ion and neutral velocity structures are highly similar and exhibit anisotropy along the local magnetic field, with a stronger anisotropy toward a larger  (see Fig. A1), similar to the anisotropy of MHD turbulence seen in a single fluid (Cho & Vishniac 2000;Cho et al. 2002;Xu et al. 2019c).The ion and neutral velocity spectra for the cases of  d = 1 × 10 5 and 1 × 10 4 follow approximately the Kolmogorov scaling with a spectral slope of −5/3, while for  d = 1 × 10 3 , the spectra become a bit steeper, with a slope of ≈ −1.9.Appendix Fig. A1 further decomposes the velocity fields within different  ranges and demonstrates the similarity in the velocity structures of ions and neutrals, irrespective of the length range.
Transition from strongly to weakly coupled regime: Fig. 2 presents the velocity distribution slices and turbulent kinetic energy spectra for three other setups with  d = 100 and 25.We find that the spectra of ions and neutrals are different, and the spectrum of ions with a slope of approximately −3.2 is steeper than that of neutrals, indicative of more severe damping of the turbulent cascade in ions.The spectrum of neutrals is also steeper than the Kolmogorov one.It suggests that the neutral-ion decoupling does not happen sharply at a particular scale, but gradually over a range of scales.Compared to the strongly-coupled case, the velocity distributions of both neutrals and ions show a clear deficiency of small-scale structures, and the anisotropy is less apparent.We can also see that the neutral-ion decoupling does not happen at the ambipolar diffusion wavenumber  dec, ∥ .Instead, only at  dec,⊥ , the spectra of ions and neutrals start to diverge.In Fig. A2, we decompose the ion and neutral velocity fields within different ranges of  for the case with  d = 25 and find that differences in velocity structures appear starting from large  values.
We note that previous studies with a low ionization fraction of 1 × 10 −4 suggest a Kolmogorov spectrum of neutrals no matter whether they are coupled or decoupled from ions (Meyer et al. 2014).In our case, we have a higher ionization fraction of 0.1 and lower  d .Our result suggests that the reduced  d may cause enhanced frictional damping and thus steepening of the spectra of neutrals and ions when they are coupled at  <  dec,⊥ .At  >  dec,⊥ , with relatively higher ion inertia neutrals are not fully decoupled from ions.Consequently, neutrals cannot develop a completely independent hydrodynamic cascade and their spectrum remains steep.This is, however, constrained by the limited internal range in our current numerical simulations.We expect the neutrals spectrum would become shallower at a sufficiently large wavenumber, in which neutrals are fully decoupled from ions.We note that our theoretically calculated decoupling and damping wavenumbers are based on the Kolmogorov scaling and scale-dependent anisotropy of Alfvénic turbulence.For strongly damped MHD turbulence with a steep spectrum and insignificant anisotropy, the theoretical estimates have a large uncertainty.Additional uncertainty comes from the fluctuations in the local ionization fraction in compressible MHD turbulence, and thus the decoupling of neutrals from ions does not happen on a singlelength scale.We further discuss this point in § 4.1.2.
Transition from weakly coupled to decoupled regime: At  d = 5, neutrals are decoupled from ions on the turbulence injection scale, while ions are still globally coupled to neutrals up to  ∼ 5 − 10.The velocity distributions of ions and neutrals exhibit differences, with neutrals displaying more isotropic velocity structures.In this regime, neutrals develop independent hydrodynamic turbulent cascades with a Kolmogorov slope, while ions undergo frequent collisions with neutrals, effectively damping the turbulence in ions.Therefore, ions exhibit a steep spectrum with a slope of approximately -2.6.The slope is related to the fraction of energy transferred to neutrals.Here we also see the ion spectrum exhibits a higher amplitude, most likely due to the artifact of driving turbulence.Our initial correlation timescale of the driving force is set to be equal to the crossing time of the Alfvén speed in the neutral and ion well-coupled cases, i.e. using the Alfvén speed calculated from the total density   +   .When neutrals decouple from ions, neutrals do not develop the Alfvén wave.The Alfvén speed then becomes larger due to smaller ion density.The correlation time fixed in the simulation is therefore too large for ions, so ions' turbulence cascades slower and gets higher velocity power.
On the other hand, in Fig. D1, we calculated the kinetic energy of ions and neutrals, the magnetic fluctuation energy, and the energy exchanged by their drag interaction.We found that the energy exchange in ions and neutrals are minimum in both  d = 10 5 and 5 cases.The velocity spectra are identical for the former, but the ion velocity spectrum has a higher amplitude for the latter.It suggests that the higher ion velocity power is not supplied by energy exchange with neutrals, but caused by the unequal correlation time of driving discussed above.

Fluctuations in local ionization fraction, Alfvén speed, and decoupling scales
The neutral-ion decoupling scale, as discussed in § 3, depends on the  ni =  d   =  d   (  +   ) (also   for the parallel decoupling scale).Due to variations in density and magnetic fields in compressible MHD turbulence, these two quantities can exhibit significant fluctuations, resulting in local variations of the decoupling scales instead of a value.
To further investigate the variation of the local decoupling scale, we present histograms of the local ionization fraction in Fig. 3 with corresponding 2D slices shown in Fig. C1.The histogram of the  d = 1 × 10 5 case is very narrow, with the ionization fraction concentrated around 0.1.However, as  d decreases, the ionization fraction starts to spread to both higher and lower values, indicating more significant local variations.For the other five cases with smaller  d , we observe that the ionization fraction varies from approximately 0 to 0.3, while the global mean value of approximately 0.1 remains the same.These variations are due to fluctuations in ion and neutral  densities.We expect that in supersonic turbulence with   much larger than unity, where density fluctuations are more significant, the variation of ionization fraction may further increase.
In addition to the ionization fraction, we also investigate the local Alfvén speed fluctuations, shown in Fig. 4 with corresponding 2D slices in Fig. C2.The Alfvén speed fluctuations come from the variation of magnetic field strength and total density   +   .Unlike the ionization fraction, the case of  d = 1 × 10 5 exhibits the widest histogram indicating significant variation of Alfvén speed.The histograms, however, become narrower for the other five cases with smaller  d .Typically, we see the maximum value of   /⟨  ⟩ reaches ∼ 2 and minimum values are either ∼ 0 (for  d = 1 × 10 4 and 1 × 10 3 ) or ∼ 0.5 (for  d = 100, 25, and 5).
The Alfvén speed and ionization fraction fluctuations result in local variations in the values of  dec, ∥ and  dec,⊥ .The distributions of their theoretically expected values calculated by using the local   and   are shown in Fig. 5, which highlights the significant fluctuations that can occur.In the case of  d = 1 × 10 5 and 1 × 10 4 , the minimum values of  dec, ∥ and  dec,⊥ are larger than the numerical dissipation wavenumber, suggesting that neutrals and ions remain locally well-coupled.Otherwise, if the local decoupling wavenumbers are smaller than the numerical dissipation wavenumber, neutrals are locally decoupled from ions.As seen in the case of  d = 1 × 10 3 , the local  dec, ∥ can vary from ≈ 1 to ≈ 1 × 10 3 , indicating the existence of local decoupling.The local  dec,⊥ can even reach a larger value of ≈ 8 × 10 3 .Neutrals can fully decouple from ions only at wavenumbers larger than the maximum  dec,⊥ .
For  d = 100, although the expected global mean decoupling scales are  dec, ∥ ≈ 10 and  dec,⊥ ≈ 32, the local values of  dec, ∥ and  dec,⊥ are also widely distributed from ≈ 1 to ≈ 50.When  d = 25, the range of  dec, ∥ and  dec,⊥ is ≈ 1 to ≈ 10, and the damping of MHD turbulent cascade in neutrals due to local coupling with ions is more noticeable than the  d = 100 case.In this case, the neutral spectrum follows the steep spectrum of ions up to  ∼ 10, which is also the maximum value seen in the histogram of  dec,⊥ .Finally, in the case of  d = 5,  dec, ∥ and  dec,⊥ do not exceed 1.5 at their maximum, indicating that neutrals are fully decoupled from ions basically at all scales and develop a hydrodynamic turbulent cascade independently.
These results suggest that neutral-ion decoupling does not occur on a single-length scale, but rather over an extended range of scales.

Density statistics
The local variation of the ionization fraction   is important to understand the neutral-ion decoupling and damping of MHD turbulence.  is directly related to the density fluctuations in ions and neutrals.In this section, we investigate the density statistics of ions and neutrals.

2D density distribution and density spectrum
Figs. 6 and 7 present 2D density slices (taken at  = 240 cell, perpendicular to the mean magnetic field) and density spectra for ions and neutrals.When  d = 1×10 5 , we observe that the density distributions of ions and neutrals are nearly identical, with the structure regulated by turbulence anisotropy.Similar filamentary density structures are also seen in single-fluid MHD simulations (Xu et al. 2019c).The spectra are a bit shallow, but they generally follow the Kolmogorov scaling, similar to their velocity spectra.However, when  d = 1×10 4 , the ion density distribution becomes different from that of neutrals.Ion density structures exhibit more apparent striations, while such small-scale structures are not seen in neutral density distribution.Correspondingly, the spectrum of the ion density becomes shallower (slope ≈ −1.1), while that of the neutral density starts to become steeper at large .These phenomena are more pronounced in the case of  d = 1 × 10 3 , where the slope of the ion density spectrum is ≈ −1.3.
We see that despite the similar velocity structures seen in neutrals and ions in the strongly coupled regime, their density structures can differ significantly.The velocity field is likely to be dominated by incompressible Alfvénic turbulence, while density fluctuations are mainly induced by compressible turbulent motions.This can be seen from the difference in the velocity and density spectra of ions.Although neutrals are strongly coupled to the Alfvénic turbulent motions, they may be poorly coupled to the compressible MHD turbulent motions and thus do not exhibit the small-scale density structures created by the compressible MHD turbulent motions.Furthermore, when  d = 100 and 25, the damping of MHD turbulence occurs (see Fig. 2).The decoupling of neutrals from the Alfvénic turbulent motions also contributes to the difference in the density structures of neutrals and ions.The density distribution in neutrals appears isotropic.We observe that the anisotropic filamentary structures in ions become less apparent, which is due to the severe damping, while sharp-density jumps gradually appear on large scales.These sharp jumps are most significant in the neutral density.The spectra of both ions and neutrals are steep for  d = 100.Together with  d = 1 × 10 3 , these three cases are complicated because of the large variation of the local decoupling scale (see Fig. 5).However, for  d = 5, the full decoupling of neutrals from ions is achieved, and only turbulence in ions is damped (see Fig. 2).In this case, we clearly see that small-scale density structures arise in neutrals and the anisotropic filamentary structures in ions vanish, and the spectra are shallower (slope ≈ −2.6 for ions and ≈ −2.1 for neutrals) than those at  d = 100 (slope ≈ −3.4 for ions and ≈ −2.7 for neutrals) and  d = 25 (slope ≈ −3.9 for ions and ≈ −2.3 for neutrals).

The probability distribution function of the logarithmic mass density
We present the probability distribution function (n-PDF) of the logarithmic mass densities of ions and neutrals, normalized by their respective mean values, as shown in Fig. 8.The n-PDF is a widely used tool for studying density statistics in single-fluid MHD turbulence (Price et al. 2011;Burkhart 2018), as it directly reveals the significance of density fluctuations.In general, the minimum and maximum values of the neutrals' n-PDF are approximately -1 and 1, respectively.These values vary a bit at large  d = 1×10 5 and 1×10 4 .As in the case of single-fluid turbulence, we expect the width of the n-PDF to be correlated with   (Padoan & Nordlund 2002).A high sonic Mach number, especially greater than unity, is typically asso-ciated with shocks, which lead to high-density contrasts and a more dispersed n-PDF.This behavior is commonly observed in studies of single-fluid MHD turbulence (Price et al. 2011;Burkhart 2018).
When  d = 1 × 10 5 , the ions' n-PDF closely resembles that of neutrals.However, the ions' n-PDFs become more dispersed with smaller  d .While the maximum value of the ions' n-PDFs remains stable at log(/⟨⟩) ≈ 1.25 − 2.0, the minimum value reaches log(/⟨⟩) ≈ −3.0 for  d = 1 × 10 3 and 100.This suggests that the ion density exhibits significant local fluctuations.The ions' n-PDFs narrow again for  d = 25 and 5, with log(/⟨⟩) ≈ −1.25 at the minimum.

Magnetic field statistics
Fig. 9 displays 2D slices of total magnetic field strength taken at  = 240 cell perpendicular to the mean magnetic field direction, as well as magnetic energy spectra calculated for the full cube.In the neutral-ion locally well-coupled state, where  d = 1×10 5 and 1×10 4 , the magnetic field fluctuations elongate anisotropically along the magnetic field direction, akin to the velocity and density structures shown in Figs. 1 and 6.The spectra exhibit Kolmogorov scaling overall.However, when  d = 1×10 3 , the magnetic field structures are  The damping of MHD turbulence 11 less filamentary, and the spectrum becomes steeper with a slope of ≈ −2.8 than that of the velocity spectra.This steepening of the magnetic energy spectrum when the velocity spectra of neutrals and ions are similar has been observed also in Meyer et al. (2014).It may be attributed to the effect of local neutral-ion decoupling.Alternatively, the fast modes in MHD turbulence may get damped at  smaller than the damping scale of Aflvén modes (Xu et al. 2015(Xu et al. , 2016;;Xu & Lazarian 2017).The damping of the magnetic fluctuations generated by fast modes may result in a steeper magnetic energy spectrum than the kinetic energy spectrum that is dominated by the Alfvén modes.
In the weakly coupled regime with  d = 100 and 25, small-scale magnetic field structures are less prominent, and the spectra become even steeper, indicating that the magnetic field energy becomes concentrated on larger scales.This is naturally expected due to the severe neutral-ion collisional damping.However, when the neutral-ion decoupling occurs at the injection scale (i.e.,  d = 5), the situation changes.The spectrum becomes shallower compared to the cases of  d = 1 × 10 3 , 100, and 25 (slope ≈ −2.8, −3.6, and −4.0, respectively).The slope is close to −2.23.It suggests a weak damping effect, as seen in Fig. 2. Overall, we see that the magnetic energy spectrum has a similar shape as the turbulent kinetic energy spectrum in ions.

Comparison with earlier studies
The two-fluid (neutral and ion) simulation requires a very short time step to stably accommodate the fastest wave speed in the problem (Meyer et al. 2014).This is computationally expensive since the low ionization fraction in the ISM (typically 10 −4 − 10 −7 for cold molecular clouds, see Tielens 2005;Draine 2011) results in an extremely large Alfvén speed.The practical time step is even smaller because the Alfvén speed in a simulation has its own variations (see Fig. C2).The stable time step is therefore determined by the largest Alfvén speed.Consequently, this computational challenge limits the numerical study of two-fluid MHD turbulence.In the study of two-fluid MHD turbulence, some earlier research has focused on the damping of MHD turbulence using simulations generated by the RIEMANN code (Tilley & Balsara 2010;Meyer et al. 2014;Burkhart et al. 2015).These studies have primarily concentrated on supersonic (  > 1) In our study, we used the Athenak code to conduct trans-sonic (  ≈ 1) two-fluid simulations with a moderately low ionization fraction of 0.1.To further examine the effect of ionization fraction, a comparison with a lower ionization fraction of 0.01 is presented in Fig. B1.Compared to earlier studies using a fixed  d , this approach of using reduced and varying  d enables us to study MHD turbulence in different coupling regimes with far fewer computational resources.Compared to earlier studies, we report newly discovered properties of two-fluid MHD turbulence.We quantitatively compared our numerical measurements with previous theoretical predictions on neutral-ion decoupling and damping scales (Xu et al. 2015(Xu et al. , 2016)), and found their large variations due to the large density fluctuations of ions and neutrals.The computational tools' application to other physical environments, like the very local ISM, will be explored in our future work.

Cosmic ray transport
The damping of MHD turbulence is crucial for understanding the transport of cosmic rays (CRs) in the multi-phase ISM.The damping of MHD turbulence should be taken into account for both resonant scattering (Xu & Yan 2013;Hu et al. 2021b) and non-resonant mirroring (Lazarian & Xu 2021).The severe damping of MHD turbulence in a weakly ionized medium can significantly affect the efficiency of scattering and the spatial confinement of CRs (Xu et al. 2016).

Velocity gradient
Determining the scale at which neutrals and ions become decoupled is a challenging task in observations.In particular, it is non-trivial to obtain the velocity spectra of ions and neutral.However, this can be achieved by the Velocity Gradient Technique (VGT; Lazarian & Yuen 2018;Hu et al. 2018), which is a new approach to tracking magnetic fields using spectroscopic data.VGT is based on the anisotropy of MHD turbulence, where turbulent eddies align themselves along the magnetic fields.Velocity gradient serves as a detector of the anisotropy and, therefore, can reveal the magnetic field direction.
The study shows that this anisotropy is absent when neutrals and ions become decoupled.We can expect that at a length scale larger than the decoupling scale, neutral turbulence and ion turbulence act as a single fluid and exhibit anisotropy, with velocity gradients of both species oriented in the same direction.At smaller length scales where neutrals decouple from ions, their velocity fields change (see Fig. 2) so that the relative orientation of their velocity gradients changes.Therefore, comparing the directions of the (ions and neutrals) velocity gradients at different length scales can independently reveal the neutral-ion decoupling scale (i.e., perpendicular ambipolar diffusion scale).This approach could provide unique constraints on the important perpendicular ambipolar diffusion scale in observation for a better understanding of star formation (Mestel & Spitzer 1956;Nakano & Tademaru 1972).

SUMMARY
Magnetized turbulence is ubiquitous in the partially ionized ISM.The interaction between neutral and ionized species can modify the properties of MHD turbulence and cause neutral-ion collisional damping.On the basis of the two-fluid MHD turbulence simulations generated from the AthenaK code, we numerically studied the statistical properties of velocity, density, and magnetic field in different regimes of ion-neutral coupling.Our main findings are: (i) Our results demonstrate that in the (neutral-ion) strongly coupled regime, velocity statistics in the two-fluid simulations can resemble those in single-fluid MHD turbulence.
(ii) In the weakly coupled regime, we observe that damping of turbulence can occur in both, resulting in their steep kinetic energy spectra compared to the Kolmogorov spectrum, while the damping of turbulence in ions is more severe.We find that due to the large density fluctuations in ions and neutrals, the ionization fraction has a large spatial variation, which causes a significant local variation of the neutral-ion decoupling scale.As a result, both neutral-ion decoupling and damping of MHD turbulence happen over a range of length scales.
(iii) In the transition regime from weakly coupled to decoupled regime ( d = 5), the damping of MHD turbulence takes place in ions showing a steep kinetic energy spectrum.Neutrals develop an independent hydrodynamic turbulent cascade and the corresponding kinetic energy follows the Kolmogorov scaling.
(iv) We find that in the strongly coupled regime with similar velocity structures in ions and neutrals, their density structures can exhibit The magnetic field maps are normalized by the mean value.The view direction is perpendicular to the mean magnetic field, which is along the vertical direction.2 nd and 4 th : magnetic energy spectrum.The spectra are averaged over several snapshots within one eddy turnover time.The shadowed areas represent the variations.significant differences.The small-scale enhanced density fluctuations seen in ions are absent in neutrals, and the density spectrum of ions is also much shallower than that of neutrals.This may be caused by the poor coupling of neutrals to the compressive MHD turbulent motions.
(v) We show that the probability distribution function (n-PDF) of neutral mass density is insensitive to the coupling status between ions and neutrals.The n-PDF of the ions is broader than that of the neutrals, and its width varies in different coupling regimes.
(vi) Finally, we find that the magnetic energy spectrum in general has a similar shape as the kinetic energy spectrum of ions.

Figure 1 .
Figure 1.Top and middle: 2D slices (taken at  = 240 cell) of ions' (top)  and neutrals' (middle) velocity field.The velocity maps are normalized by the mean value.The view direction is perpendicular to the mean magnetic field, which is along the vertical -direction.Bottom: Turbulent kinetic energy spectra of ions (red) and neutrals (blue).The spectra are averaged over several snapshots after turbulence reaches a statistically steady state, with the time interval equal to the largest eddy turnover time.The shadowed areas represent the variations.Magenta and green lines represent the theoretically expected neutral-ion parallel decoupling (Eq.5) and damping wavenumbers (Eq.6), respectively, for Alfvén modes of MHD turbulence.The perpendicular decoupling and damping wavenumbers are not shown because they are larger than 10 3 .

Figure 6 .
Figure 6.Top and middle: 2D slices (taken at  = 240 cell) of ion (top) and neutral's (middle) density field.The density maps are normalized by the mean value.The view direction is perpendicular to the mean magnetic field, which is along the vertical direction.Bottom: Spectra of ion (red) and neutral's (blue) density.The spectra are averaged over several snapshots with one latest eddy turnover time.The shadowed areas represent the variations.

Figure 8
Figure 8. n-PDFs of ion (red) and neutral's (blue) logarithmic mass densities normalized by their mean densities.

Figure 9 .
Figure 9. 1 st and 3 rd : 2D slices (taken at  = 240 cell) of magnetic field strength.The magnetic field maps are normalized by the mean value.The view direction is perpendicular to the mean magnetic field, which is along the vertical direction.2 nd and 4 th : magnetic energy spectrum.The spectra are averaged over several snapshots within one eddy turnover time.The shadowed areas represent the variations.

Figure A1 .Figure A2 .
Figure A1.2D slices (taken at  = 240 cell) of ion velocity (top) and neutral velocity (bottom) for the simulation of  d = 1 × 10 3 .The velocity fields are decomposed into different  ranges in Fourier space and then transformed back to real space.

Table 1 .
Setups of two-fluid simulations.  and   are the instantaneous RMS values at each snapshot that is taken. = 2(   /  ) 2 is plasma compressibility. dec,∥ and  dec,⊥ are theoretically expected parallel and perpendicular components of the neutral-ion decoupling wavenumber, respectively.The listed  d is given in numerical units.To obtain a dimensionless value, divide  d by  inj /(  inj   ), which is fixed at 10 and 100 (in numerical units) for   = 0.1 and 0.01, respectively.
dis ≈ 40−50, indicating that ions and neutrals are well-coupled over all length scales resolved in our simulations.Note  d values are given in numerical units.To obtain a dimensionless value, one can divide  d by  inj /( inj   ), which is approximately 10 for the simulations with   = 0.1.To calculate the theoretically expected  dec, ∥ and  dec,⊥ , as well as  dam, ∥ and  dam,⊥ , we adopt the mean values of density, magnetic field, and  inj at  = 1.We denote the  dec, ∥ and  dec,⊥ as averaged decoupling wavenumbers.If the averaged decoupling wavenumbers are larger than the numerical dissipation wavenumbers, then neutrals and ions are on average well-coupled.