© The Authors 2025

1 Introduction

Recent high-resolution observations of young stellar objects with the Atacama Large Millimeter/submillimeter Array (ALMA) at (sub)millimeter wavelengths reveal several substructures in their circumstellar disks, such as rings, spirals, or gaps (e.g., ALMA Partnership 2015; Andrews et al. 2018). These substructures are either attributed to still evolving and forming but yet unseen giant planets or brown dwarfs (Brown et al. 2009; Andrews et al. 2011; Biller et al. 2012), or to magneto-hydrodynamical instabilities that play an important role during the accretion process.

In the inner disk regions, that is, inside 1 au, the magneto-rotational instability (Balbus & Hawley 1991; Flock et al. 2011; Turner et al. 2014) and magnetically driven winds (Bai et al. 2016; Gressel et al. 2020; Lesur 2021) determine the accretion flow and shape of the disk structure (Flock et al. 2017a; Lesur et al. 2023). In contrast, for less ionized regions in the outer disk, hydrodynamical instabilities might become more important in determining the gas and dust flows. Most notably, the vertical shear instability (VSI, Urpin & Brandenburg 1998; Urpin 2003; Arlt & Urpin 2004; Nelson et al. 2013; Stoll & Kley 2014; Lin & Youdin 2015; Flock et al. 2017b; Schäfer et al. 2020; Flores-Rivera et al. 2020; Barraza-Alfaro et al. 2021), which is related to the Goldreich–Schubert–Fricke instability (Goldreich & Schubert 1967; Fricke 1968), is a promising mechanism that causes turbulence in the regions outside 10 au of protoplanetary disks. This instability arises from fluctuations in the orbital motion in the vertical direction, that is, along the rotational axis of the disk.

Additionally, it has been shown that the VSI can cause the formation of vortices, which act as particle traps and are therefore essential in the context of planet formation (Richard et al. 2016; Latter & Papaloizou 2018; Manger & Klahr 2018). In this work, we use simulation results from high-resolution 3D radiation hydrodynamical simulations of the outer regions of protoplanetary disks presented in Flock et al. (2020). This work also includes the individual motion of larger grains as a result of the gas drag.

Once interstellar micrometer-sized grains reach higher densities in the disk, they start to grow and settle to the midplane of the disk (Beckwith et al. 2000). Micrometer-sized grains are still well coupled to the gas motion. However, once they grow to sizes of hundreds of micrometers, their motion starts to decouple from that of the gas phase (Weidenschilling 1977). During this phase of dust evolution, a large amount of solid material, in the form of millimeter- and centimeter-sized pebbles, is being rearranged in the young protoplanetary disk, mainly by the fast inward radial drift and the vertical migration toward the disk midplane (Birnstiel et al. 2010; Andrews 2020; Birnstiel 2024).

These irregularly shaped dust pebbles typically align on large scales along a preferential direction relative to a given axis due to torques acting upon them. Such torques can arise from radiation, magnetic fields, or mechanical gas drag. The most promising alignment mechanism is radiative torque (RAT) alignment (e.g., Dolginov & Mitrofanov 1976; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007a,b; Hoang & Lazarian 2008). In this scenario, radiation spins up the dust grains, and the rotation axis subsequently aligns with the magnetic field orientation due to paramagnetic dissipation, a mechanism that has been extensively studied in the context of circumstellar disks (Cho & Lazarian 2007; Reissl et al. 2016; Tazaki et al. 2017). As a result, the net thermally re-emitted radiation from the aligned, elongated grains is expected to be partially polarized, with the observable polarization vectors tracing the magnetic field orientation (see, e.g., Andersson 2015, for a review).

However, the question of which mechanisms, other than RATs, can drive coherent grain alignment remains an active area of research (Dolginov & Mitrofanov 1976; Purcell 1979; Andersson 2015; Hoang et al. 2018). A mechanism based on supersonic gas-dust drift was proposed more than half a century ago (Gold 1952a,b). Such a gas-dust drift would mechanically align spheroidal grains by minimizing their geometrical cross section. In this scenario, the resulting dust polarization is related to the direction of the drift velocity and does not trace the magnetic field. A description of mechanical alignment based on a toy model was provided by Lazarian & Hoang (2007a), Das & Weingartner (2016), and Reissl et al. (2023), suggesting that irregularly shaped dust grains may also attain a long-term stable alignment configuration even under subsonic drift conditions. In the case of mechanical alignment, the helicity of an oblate grain can be either right-handed or left-handed (Kataoka et al. 2019). In both cases, the resulting polarization is oriented perpendicular to the direction of the gas velocity.

Consequently, measuring the degree of polarization of dust thermal emission is a potentially powerful tool for constraining the physical processes related to gas flow in protoplanetary disks. Therefore, in this study, we do not attempt to reproduce previous efforts to model dust polarization in disks based on RATs and magnetic field alignment (e.g., Cho & Lazarian 2007; Tazaki et al. 2017), but instead focus on the processes that enable the mechanical alignment of dust pebbles. Specifically, the goal of this study is to analyze the impact of mechanical grain alignment caused by the drift motion of these pebbles in the protoplanetary disk on the resulting observable linear polarization.

We used simulation snapshots from 3D radiation hydrodynamical simulations, including the motion of larger dust pebbles (Flock et al. 2020). Subsequently, we constrained the characteristics of wavelength-dependent polarization maps resulting from the polarized emission of aligned grains. Here, we consider synthetic observations in the millimeter wavelength range, that is, from 0.87 mm to 10 mm, to trace the thermal emission of large millimeter-sized dust grains.

This paper is organized as follows. In Sect. 2, we introduce and describe our model setup and numerical procedure. Subsequently, we present and discuss synthetic observations generated using radiative transfer simulations in Sect. 3. In Sect. 4, we discuss several assumptions and limitations in our grain alignment model, and finally we summarize our study in Sect. 5.

2 Model setup

To constrain the orientation of drifting grains under the effect of mechanical alignment, we applied 3D radiation hydrodynamical simulations (Flock et al. 2017b, 2020) using the PLUTO code (Mignone et al. 2007, 2010, 2012). Here, simulation snapshots provide the density of the gas and dust, as well as the velocity fields of the gas v_g and the dust component v_d in each cell. We then define the drift velocity or headwind of the dust grains as $${\text{SO}}_4^{2-}$$(1)

In this study, we used the data of the simulation snapshot after 750 orbits of evolution. After 650 orbits of evolution, particles of 0.1 mm and 1 mm in size were inserted at the midplane. We let the particles evolve for a further 100 orbits to allow them to drift and mix into a new equilibrium state. An overview of the setup parameter is described in Sect. 2.4 and the corresponding Table 1. For more information on the solved equations and the numerical setup, we refer to Sect. 2 in Flock et al. (2020).

2.1 Ballistic grain aggregation

Modeling the mechanical alignment of dust is closely related to the shape and complex surface topology of individual grains. Observations and modeling of (sub)millimeter polarization suggest the presence of compact grains in the disk midplane (Stephens et al. 2020; Brunngräber & Wolf 2021), while scattering properties derived from mid-infrared observations suggest more elongated porous grains (Ginski et al. 2023; Tanaka et al. 2023). The porous grain aggregates are likely to undergo some compression while settling toward the midplane (Tanaka et al. 2023; Ueda et al. 2024).

In this paper, we utilize moderately compact dust grains throughout the entire disk to estimate the parameters of mechanical alignment. Dust grain aggregates are created using ballistic aggregation and migration (BAM) of primary particles (monomers). Here, we follow the algorithm suggested in Shen et al. (2008), with modifications presented in Reissl et al. (2024), where monomers on random trajectories are successively added to the aggregate. For this study, we allow for exactly one migration (BAM1) to establish additional connections, controlling the compactness of the aggregates. The monomer radii are sampled from a log-normal distribution between 10 nm and 100 nm with a mean of 20 nm. The sampling from the log-normal distribution is biased to guarantee an exact pre-defined effective radius a_eff = (3 ∑ V_mon/(4 π))^1/3, where the total volume of the resulting aggregate is the sum of the volumes of the individual monomers V_mon.

We pre-calculate aggregates with exact effective radii between 0.1 μm and 1.4 μm. Calculating aggregates with larger effective radii a_eff is challenging due to increasing complexity within the BAM algorithm and resulting computational time constraints. For each a_eff, we use 30 random seeds to create a representative sample of BAM1 aggregates. An exemplary BAM1 dust aggregate with a_eff = 1.2 μm is depicted in Fig. 1.

Fig. 1 Exemplary BAM1 dust aggregate with an effective radius of aeff = 1.2 μm. The aggregate rotates with an angular velocity, ω, that corresponds to its maximum moment of inertia, Im. The rotation and precession of the dust aggregate result from gas-dust interactions, where the gas and dust components move with a relative velocity, vdrift, through the disk. The angle, Θ, is defined as the angle between ω and vdrift.

2.2 Mechanical grain alignment

A gas-dust drift has long been suggested as a cause of mechanical alignment (Gold 1952a,b). However, attempts to constrain the parameters of mechanical alignment lack predictive capability for a larger grain ensemble or do not sufficiently take the complex surface morphology of dust aggregates into account (see, e.g., Lazarian 1994; Lazarian & Hoang 2007b; Das & Weingartner 2016; Hoang et al. 2018).

In this paper, we apply the mechanical alignment of dust (MAD), a numerical framework for calculating alignment parameters using detailed Monte Carlo dust-gas interaction simulations. The MAD framework considers the mechanical torques (MET) arising from scattering, sticking, and desorption of gas on the surface of fractal dust aggregates (see Reissl et al. 2023, for further details). However, it is not a priory clear how the mechanical alignment of BAM1 aggregates would scale with an increasing grain size of a_eff. Hence, we explore the MAD for grain sizes of a_eff = 0.1 μm to 1.4 μm in order to derive a parametrized expression for the mechanical alignment.

In detail, the MAD is used to calculate the METs for each of the individual pre-calculated BAM1 aggregates introduced in Sect. 2.1 for representative values of the gas temperature, T_g, gas number density, n_g, dust temperature, T_d, and the differential velocities of the gas component, v_g, and the dust component, v_d, separately. The MAD simulation results are analyzed by utilizing the dimensionless quantity v_drift/v_th, that is, the absolute drift velocity, v_drift = |v_drift| normalized by the thermal velocity of the gas $$v_{\mathrm{t}\mathrm{h}}=\sqrt{\frac{2k_{\mathrm{B}}{T}_{\mathrm{g}}}{\mu m_{\mathrm{H}}}}$$(2)

Here, we assume that the gas temperature, T_g, is equal to the average dust grain temperature and an average molecular weight of μ = 2 assuming molecular hydrogen.

When modeling the rotational dynamics, the dust aggregate cannot be assumed to be a rigid body. Internal dissipative processes (see, e.g., Purcell 1979; Lazarian & Roberge 1997; Lazarian & Draine 1999; Lazarian & Efroimsky 1999) make, in the long run, the axis associated with the maximal moment of inertia, I_m, the axis of rotation. For simplicity, we assume the dissipation timescale of the dust aggregates to be small compared to the dynamical timescales of the disk. Subsequently, the gas-dust interactions spin up a dust aggregate to an angular velocity ω. We note that the METs acting on the dust aggregate would not align ω parallel to the direction of the drift v_drift but would rather cause a precession of ω around v_drift with an opening angle of Θ. A sketch of the alignment dynamics is provided in Fig. 1.

For each set of input parameters {n_g, T_g, T_d, v_drift/v_th} we evaluate the dynamical evolution of the angular velocity ω and the opening angle Θ via $$\frac{\mathrm{d}\omega }{\mathrm{d}t}=\frac{1}{I_{\mathrm{m}}}{\mathbf{\Gamma }}_{\mathrm{M}\mathrm{E}\mathrm{T}}(\mathrm{\Theta })-\frac{\omega }{{\tau }_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}}$$(3)

for each BAM1 aggregate individually. Here, Γ_MET(Θ) is the angle-dependent net MET, and τ_drag is the characteristic timescale of the rotational drag, combining the effects of gas drag and IR drag due to photon emission (Draine & Lazarian 1998). For the calculation of τ_drag of BAM1 aggregates, we refer the reader to Reissl et al. (2023) and Jäger et al. (2024), respectively. A solution to Eq. (3), along with the average grain alignment behavior as a function of grain size and the corresponding timescales, is presented in Appendix A.

Furthermore, we trace the trajectory for each of the BAM1 aggregates with the same effective radius, a_eff, in the ω−Θ phase space defined by Eq. (3) to determine attractor points that are the long-term points of stability of the angular momentum, ω, and the opening angle, Θ. As a lower limit for grain alignment, we take the thermal angular momentum ω_th ≈(k_B T_d/I_m)^1/2. Aggregates rotating below ω<3 ω_th may still have attractors, but are regularly kicked out of stable alignment by random gas-dust collisions (Hoang & Lazarian 2008) and, thus cannot contribute to the net dust polarization.

An upper limit arises from the fact that rapidly rotating dust grains may become disrupted by centrifugal forces, given a sufficiently high angular velocity ω_disr (see Hoang et al. 2019, for further details). Recently, a study by Reissl et al. (2024) based on N-body simulations of rotating aggregates provided a best-fit parameterization of ω_disr based on basic material properties. We utilize this parameterization to estimate the limit for rotational disruption, ω_disr, for each individual BAM1 grain, assuming typical material parameters of composite silicate grains as presented in Reissl et al. (2024). Using the exact grain composition introduced in Sect. 2.1 would require modeling dust aggregates where the water ice is not uniformly distributed within the aggregate but frozen to the aggregate’s surface. Such modeling with unevenly distributed materials goes beyond the scope of this paper. However, given the variations in the N-body simulation results of Reissl et al. (2024) for different grain compositions, we estimate that the magnitude of ω_disr may be off by a factor of two at most.

Finally, we introduce the ratio, ζ. This ratio quantifies the number of BAM1 aggregates that have an attractor point within 3 ω_th < ω < ω_disr relative to the total number of pre-calculated aggregates per effective radius, a_eff, to quantify the dust polarization of an ensemble of dust grains (see Reissl et al. 2023, for further details). The maximum possible dust polarization is determined by the refractive index and the shape of the dust grain. However, due to the precession of the dust grains, the maximal polarization becomes further reduced. This reduction is usually quantified by the Rayleigh reduction factor (see, e.g., Greenberg 1968; Lee & Draine 1985; Roberge & Lazarian 1999) $$\mathcal{R}=\frac{3\zeta }{2}(⟨{\cos }^2\mathrm{\Theta }⟩-\frac{1}{3}),$$(4)

where 〈cos² Θ〉 is the ensemble average of all the angles Θ of the individual BAM1 aggregates at their attractor points. We emphasize that we use a definition of the Rayleigh reduction factor modified by the ratio ζ. Exemplary distributions of Θ and the corresponding values of ℛ are presented in Appendix B for selected grain ensembles of similar size.

What the MAD Monte Carlo simulations in Reissl et al. (2023) reveal is that, by increasing the gas density, n_g, the MET, Γ_MET, increases, as does the gas drag. Beyond a certain critical density, n_c ≈ 10² cm⁻³, the MET and gas drag reach a balance, and the angular velocity ω of the aggregates reaches its terminal value. For the range of gas densities n_g ≈ 10⁷ cm⁻³ to 10¹¹ cm⁻³ present in our simulation data, each individual dust grain is already beyond this critical limit n_c. Therefore, the grain alignment process is considered to be independent of n_g for this particular PLUTO simulation.

The same applies to the gas temperature, T_g, and the dust temperature, T_d. An increase in T_g raises the impact velocity (see Eq. (2)) and subsequently increases Γ_met.. With an increase in T_d, IR photons carry away a fraction of the rotational energy, slowing down the grain rotation. However, since both temperatures, T_g and T_d, are below 100 K (compare Fig. D.1), their impact on the overall grain rotation dynamics is minimal, making the ratio v_drift/v_th the dominant parameter for modeling dust rotation and subsequent polarization.

In Fig. 2, we present the average angular momentum ω as a function of the ratio v_drift/v_th under typical disk conditions: n_g = 10⁸ cm⁻³, T_g = 20 K, and T_d = 20 K, with effective grain radii a_eff = 0.1 μm to 1.4 μm. As drift velocity, v_drift, increases, the aggregates experience collisions predominantly from the direction of v_drift rather than from random directions. This results in a systematic mechanical torque Γ_MET with a steady orientation over time. Additionally, an increase in v_drift raises Γ_MET due to a higher gas-dust collision rate. The resulting average angular momentum ω for aggregates with a_eff = 0.1 μm surpasses the threshold for stable alignment, 3 ω_th, at v_drift ≈ 0.03 v_th, and such aggregates become rotationally disrupted on average at v_drift ≈ v_th.

We emphasize that individual grains may still have stable alignment attractors within the range 3 ω_th < ω < ω_disr even though the average angular velocity, ω, does not meet this criterion. Thus, the Rayleigh reduction factor ℛ ≈ 0 to 0.75 for the ensemble of the smallest grains within the range 0.01 ≲ v_drift/v_th ≲ 3. For the largest aggregates, with a_eff = 1.4 μ m, the thermal angular momentum ω_th is lower due to their larger moment of inertia, I_m, while the limit ω_th remains in the same order of magnitude for all aggregates (see Reissl et al. 2024).

As reported in Reissl et al. (2023), larger grains are more efficiently spun up due to their increased surface area, whereas for grains with rounder shapes, the net mechanical torque cancels out over time. For the specific BAM1 aggregates used in this study, smaller grains are more elongated, while the shapes of larger grains become increasingly spherical (Reissl et al. 2024). Considering both factors – namely, the increase in surface area countered by the reduced net torque due to a rounder shape – the average behavior of the BAM1 aggregates becomes more uniform with increasing radius a_eff. What we report here is that the mechanical alignment behavior converges for grain sizes a_eff ≥ 1.0 μm. Consequently, our approach to model the mechanical alignment for individual grain sizes of a_eff = 0.1 μm to 1.4 μm allows us to extrapolate the alignment behavior up to a grain size of 100 μm and 1 mm, respectively (see Table 1). As shown in Fig. 2, this results in a Rayleigh reduction factor that can be approximated by $$\mathcal{R}=\{\begin{array}{ll}0.29& \text{\hspace{0.17em}if\hspace{0.17em}}1.78\times {10}^{-2}\le v_{\text{drift\hspace{0.17em}}}/v_{\text{th\hspace{0.17em}}}\le 10\\ 0& \text{\hspace{0.17em}otherwise\hspace{0.17em}}\end{array}$$(5)

for larger grains. For simplicity, we used this approximation of ℛ for all dust grains with a_eff ≥ 1 μm in the subsequent RT simulations. Here, the velocity fields of the gas and each dust grain in the protoplanetary disk are provided by the simulation snapshot (see Appendix C). Consequently, the direction of alignment is determined by the direction of the drift velocity v_drift, where helical grains align with their short axis parallel to the drift velocity vector.

Fig. 2 Left panel: Angular velocity, ω, as a function of the ratio vdrift/vth. The effective grain radii, aeff = 0.1 μm to 1.4 μ m, are color-coded. Solid lines represent the average values of ω, while the shaded areas indicate the range between the minimum and maximum values. Dotted lines show the lower limit, 3 × ωth, for long-term stable alignment, and dashed lines represent the upper limit for rotational disruption, ωdisr. Right panel: Same as the left panel, but for the Rayleigh reduction factor ℛ (solid lines). Dashed-dotted lines depict the parametric representation of ℛ for grains with aeff ≥ 1 μm. See Sect. 2.2 for details.

2.3 Optical properties of dust grains

To calculate the optical properties of dust grains of different sizes at various wavelengths in a reasonable amount of time, we ignore any porosity or irregularly shaped particles. In particular, calculation of the optical properties for various incident and scattering angles for the exemplary dust aggregate shown in Fig. 1 would be computationally too time-consuming and is out of the scope of the current study. We refer the reader to Kirchschlager & Wolf (2014) and Kirchschlager et al. (2019) who discussed the impact of porosity on the polarization degree, or to Kirchschlager & Bertrang (2020) who discussed the scattering of nonspherical dust grains on the polarization. Instead, we assume compact spherical grains for dust grain heating and polarization as a result of scattering and elongated oblate-shaped grains for polarization as a result of dichroic emission and absorption.

The dust grains are divided into three size bins. The first bin represents small grains with sizes of 0.1 μm to 10 μm distributed with a power law, $$N(a)\propto a^{-3.5}.$$(6)

These dust grains are fully coupled to the dynamics of the gas during hydrodynamical simulations and thus are not affected by mechanical grain alignment in this study. The other two bins represent larger grains with a log-normal grain size distribution, $$N(a)\propto \frac{1}{a}\exp (-\frac{{(\ln \left(a/a_{\mathrm{e}\mathrm{f}\mathrm{f}}\right))}^2}{2{\sigma }^2}),$$(7)

where σ=0.25 and a_eff is either 100 μm or 1 mm. These grains are less coupled to the gas movement, thus, experience a headwind from the gas and are able to align in the protoplanetary disk. We note that in the 3D hydrodynamical simulations, only single grain sizes of 100 μm and 1 mm are applied. For polarized radiative transfer, we take a size distribution since dust grains in circumstellar disks occur in different sizes and, thus, contributions of single grains are smoothed or cancel out. In Fig. 3, we present the dust mass density distribution of the different grain sizes resulting from a simulation using the PLUTO code. The wave pattern of the large dust grains is due to the large-scale upward and downward gas movement by the VSI (Flores-Rivera et al. 2020).

We apply the DSHARP dust composition (Birnstiel et al. 2018), consisting of refractory organics (Henning & Stognienko 1996) with a fraction of 40% by mass, 33% astronomical silicate (Draine 2003), 20% water ice (Warren & Brandt 2008), and 7% troilite (Henning & Stognienko 1996). The composition is mixed using the formula of Bruggeman (1935), and the resulting material density is 1.675 g cm⁻³.

For spherical particles, the optical properties of dust grains are calculated using the Mie scattering theory (Mie 1908) and the code MIEX (Wolf & Voshchinnikov 2004, 2018). Optical properties for nonspherical 100 μm and 1 mm-sized dust grains are calculated using the discrete dipole approximation (Draine & Flatau 1994) and the code DDSCAT 7.3.3 (Draine & Flatau 2000, 2008)¹. Here, we assume that the spheroidal particles are oblate with a ratio of the minor to major axes of 1/1.5. For large dust grains, that is, a ≳ λ/2, where λ is the wavelength of radiation, the polarized emission becomes unpolarized even if the grains are aligned (Cho & Lazarian 2007; Kirchschlager et al. 2019). Therefore, if a/λ ≥ 1, instead of using DDSCAT, we assume spherical particles and calculate the optical properties with MIEX. For the largest 1 mm-sized grains, we obtain an average absorption mass opacity of approximately 0.78 cm² g⁻¹ and an average scattering mass opacity of approximately 11.9 cm² g⁻¹ at a wavelength of 3 mm.

Figure 4 shows the absorption efficiency of an oblate spheroid as a function of radius of an equal volume sphere at a wavelength of 3 mm. If the grain size is small compared to the wavelength, that is, if a/λ ≲ 0.1, the ratio of the absorption efficiency is approximately 1.8. With increasing grain size, that is, 2 π a ≳ λ, the cross sections of single grain sizes show oscillations due to multiple internal reflections of radiation (van de Hulst 1981). This also causes a ratio of absorption efficiency below 1 and a change in the orientation of polarization (see Cho & Lazarian 2007; Kirchschlager et al. 2019; Guillet et al. 2020). Moreover, porosity of dust grains would reduce the intensity of oscillations in the cross sections in the Mie regime (2 π a ∼ λ) compared to compact dust particles.

Fig. 3 Cross-section view of the dust mass density in units of g cm−3 of 0.1 μm to 10 μ m-sized grains (green), 100 μ m-sized grains (purple), and 1 mm-sized grains (orange). The values are azimuthally averaged.

Fig. 4 Ratio of absorption efficiencies Q⊥/Q‖(top) and absorption efficiency (bottom) of an oblate spheroid as a function of radius of an equal volume sphere. This figure shows an axis ratio of the spheroid of 1/1.5 and a wavelength of radiation of 3 mm. The absorption efficiency Q‖is the efficiency along the minor grain axis, that is, the axis of alignment. If a/λ ≥ 1, we assume a spherical particle, thus, Q⊥ = Q‖. See Sect. 2.3 for details.

2.4 Radiative transfer

The prediction and analysis of multi-wavelength polarization measurements requires a comprehensive understanding of radiation transfer, dust properties, and underlying physical mechanisms such as grain alignment. To model synthetic observations of protoplanetary disks, we apply the 3D Monte Carlo radiative transfer code POLARIS² (Reissl et al. 2016, 2018). The model parameters are summarized in Table 1. The stellar parameters are adopted from Flock et al. (2020) and Blanco et al. (2021).

First, based on the dust density distribution, POLARIS calculates the temperature of the dust grains that are heated by the central stellar radiation source. The resulting temperature distribution is shown in the Appendix D. Subsequently, the polarized radiation from dichroic emission and absorption for aligned oblate spheroids was calculated. However, the dust grains are not perfectly aligned, causing a decrease in the observed net polarization compared to that expected for the individual grains. To account for this effect, we use the Rayleigh reduction factor defined in Eq. (5). In this work, we consider either 100 μm or 1 mm grains aligned, and we study each of their contributions to the polarization at different wavelengths. The respective other grain size, as well as the smaller grains (0.1 μm to 10 μm) have a random orientation. Finally, the polarized scattered radiation of the star and dust, that is, self-scattering, is computed by assuming spherical dust grains. In this step, the contribution of all grains to the net polarization is taken into account.

To describe the state and degree of polarization of the radiation, we use the Stokes formalism (Bohren & Huffman 1998). Here, the radiation is characterized by its total flux I, the components Q and U, which describe the linear polarization, and the component V for the circular polarization. The degree P and orientation ψ of linear polarization are defined by $$P=\frac{\sqrt{Q^2+U^2}}{I},\tan (2\psi )=\frac{U}{Q}.$$(8)

For a detailed description of the computational methods, we refer the reader to Reissl et al. (2016) and Appendix B in Reissl et al. (2019).

3 Results

Before calculating the flux and polarization maps, we have to determine whether grains are aligned in the protoplanetary disk and, in that case, the orientation of their axis of alignment. Most importantly, this is based on the drift velocity and temperature at each location in the protoplanetary disk. Subsequently, we analyze the contributions of the thermally reemitted radiation and dichroic extinction of the aligned grains as well as of self-scattered radiation to the polarization.

3.1 Orientation of aligned grains

Figure 5 shows the ratio of the drift velocity and the thermal velocity (see Eq. (5)) of 100 μm and 1 mm-sized grains in the disk midplane. Regardless of the grain size and distance to the center, most dust grains meet Eq. (5) and are therefore aligned due to METs. The minor axis of the oblate grain is then oriented parallel to the drift velocity vector v_drift. Closer to the inner rim of the disk, there is a pattern of changing velocity ratio v_drift/v_th, in particular for grains with a size of 1 mm. This can be traced back to an alternating velocity field close to the inner radius of the disk (see Appendix C) and is a numerical effect, as the region inward 25 au is affected by the boundary conditions. In contrast, an interesting feature is the change of alignment of grains trapped in a vortex. There, the velocity ratio decreases to the lower limit of Eq. (5), in contrast to the local environment.

Due to the VSI, the dominant component of the drift velocity is parallel to the rotational disk axis (see Appendix C). Thus, most of the elongated grains are oriented with their short axis perpendicular to the disk midplane, since we assume that the grains align with their short axis parallel to the drift velocity vector. As a consequence, the polarization is close to zero for a face-on disk since, in this case, the main axis of grain alignment is parallel to the line of sight, and oblate-shaped particles appear symmetric for this viewing geometry. Therefore, we consider an inclined disk for the following analysis.

3.2 Polarization of aligned grains

Figures 6 a and 6 b show the surface brightness map of the thermal reemission at wavelengths of 1.3 mm and 10 mm, respectively. The disk has an inclination of 45^∘. We assume that at wavelengths of 1.3 mm and 10 mm, only grains with a size of 100 μm and 1 μm are aligned due to METs and contribute to the polarization, respectively.

In general, the thermal reemission flux decreases with increasing distance to the center and with increasing wavelength. Due to the radially inhomogeneous distribution of dust grains, there are gaps and rings visible on the surface brightness maps. Furthermore, there is a local increase in surface brightness caused by the accumulation of dust grains due to a vortex formed at a position angle of approximately 225^∘. At a wavelength of 10 mm, the contrast in surface brightness with the surroundings increases even more and is a prominent feature there. We also refer the reader to Blanco et al. (2021) who discussed the impact of substructures produced by the VSI on synthetic photometric observations.

Figure 6 also shows the degree and orientation of polarization. Hereby, the direction of polarization is parallel to the major axis of the disk. This can be explained by the drift velocity vector, which has its largest component in the polar direction, that is, perpendicular to the disk midplane (see Sect. 3.1).

In the next step, we investigate the radial dependence of the degree of polarization. For this purpose, the polarization map of the inclined disk was deprojected and then averaged azimuthally to obtain the radial profile (see Fig. 7). This is justified since the polarization profile along the major and minor axes is identical for our model setup and the considered wavelengths.

If only 100 μm-sized grains are aligned, the degree of polarization is small throughout the disk, which are values below 1%. This is due to the decreasing absorption efficiency with decreasing dust grain size (see Fig. 4), thus, the thermal emission and polarized radiation is dominated by the larger 1 mm-sized grains. In addition, the degree of polarization increases with increasing distance to the center. This finding can be explained by the large millimeter-sized grains that are settled to the midplane and migrated inward to the inner rim of the disk (see Fig. 3).

If only 1 mm-sized grains are aligned, the degree of polarization significantly depends on the considered wavelength. In particular, for a wavelength of 10 mm, the degree of polarization increases to about 2.5%. In contrast, for a wavelength of 1.3 mm, the degree of polarization is below 0.2% throughout the disk. Moreover, the radial polarization profile decreases with increasing radius for wavelengths of 3 mm and 10 mm. This can also be explained by the large millimeter-sized grains that settled to the midplane and migrated inward to the inner rim of the disk. As a consequence, the polarization due to thermally emitted radiation decreases with increasing radius.

In addition, there are multiple maxima and minima in the polarization profiles. These are caused by gaps and rings in the surface brightness, and thus minima and maxima in the dust density distribution. Depending on the total optical depth of the structure and thus the observing wavelength, a gap or bright rings cause either an increase or decrease in the total degree of linear polarization. Most interestingly, at a wavelength of 7 mm, the degree of polarization also decreases closer to the inner rim of the disk, which is in contrast to the outlined behavior at wavelengths of 3 mm and 10 mm. Finally, at a wavelength of 10 mm and a radius of about 30 au, the degree of polarization has a local minimum. This is the result of the formed vortex, which has a lower degree of polarization compared to regions with the same distance to the center.

Fig. 5 Ratio of the drift velocity and the thermal velocity in the disk midplane of 100 μm (top) and 1 mm-sized (bottom) grains. Dust grains are aligned if 1.78 × 10−2 ≤ vdrift/vth ≤ 10 (see Eq. (5)), otherwise they are marked in red in the figure. See Sect. 3.1 for details.

3.3 Polarization due to self-scattering

Polarization arises not only from emission by aligned grains, but also from radiation emitted and scattered by dust grains before leaving the disk (self-scattering). Figures 6c and 6d shows the surface brightness map of the thermal radiation at a wavelength of 1.3 mm and 10 mm, respectively. Again, the disk has an inclination of 45^∘. Polarization is caused only by self-scattered radiation, whereas grains of all sizes contribute to the net polarization.

In general, scattered radiation contributes significantly to the observed polarization at these wavelengths, resulting in a degree of polarization of about ≲ 1.5% for all considered wavelengths. Here, the highest degree of polarization is obtained for the shortest considered wavelength of 0.87 mm. With increasing wavelength, thus decreasing optical depth, the contribution of self-scattering decreases as well. However, since most 1 mm-sized particles have already migrated inward, polarization due to self-scattering has the greatest contribution at the inner disk rim, in particular at a wavelength of 10 mm (see Fig. 6).

The orientation of the resulting polarization is usually parallel to the projected minor disk axis. This is caused by the anisotropic radiation field and the fact that scattered radiation is polarized perpendicular to the scattering plane, in particular, for grains that are smaller than the wavelength (Rayleigh limit). On the other hand, for a face-on disk, the polarization vectors would preferentially be averaged azimuthally. We refer to Kataoka et al. (2015) and Yang et al. (2016) who discuss this mechanism for thermal radiation in detail. Note that the effect of polarization reversal (see Sect. 3.5), as discussed by Yang et al. (2016) and Brunngräber & Wolf (2019), may cause a radial polarization pattern.

Figure 7 also shows the wavelength dependence of the radial profile of the degree of polarization due to self-scattering. Again, to derive this profile, the polarization map of the inclined disk was deprojected and then azimuthally averaged. We note that in contrast to the polarization of thermal emission, the polarization due to self-scattering slightly differs between the projected major and minor axes of the protoplanetary disk. In particular, the degree of self-scattering polarization is larger along the major axis of the disk compared to the minor axis of the disk (see Yang et al. 2016).

Most importantly, in Fig. 7, the degree of polarization of self-scattered is in the range or even exceeds the polarization due to emission of aligned 100 μm-sized grains. As a consequence, polarization due to self-scattering is the dominant polarization mechanism in this case. In contrast, polarized emission of aligned grains with a size of 1 mm, in particular at wavelengths of 7 mm and 10 mm, exceeds the polarization due to self-scattered radiation. Furthermore, the polarization profile decreases with increasing radius at wavelengths of 7 mm and 10 mm. Similarly to polarized emission, this is caused by large millimeter-sized grains that settled down to the midplane and migrated inward to the inner rim of the disk (see Fig. 3).

Fig. 6 Surface brightness map of the total emitted thermal radiation at a wavelength of 1.3 mm (left) and 10 mm (right) overlaid with polarization vectors at these wavelengths. In the top row (a, b), the polarization arises from thermal emission and dichroic extinction of aligned dust grains only. At wavelengths of 1.3 mm and 10 mm, only 100 μm and 1 mm-sized grains are assumed to be aligned, respectively. In contrast, in the bottom row (c, d), the polarization arises from self-scattered radiation only. The disk is inclined by 45∘. See Sects. 3.2 and 3.3 for details.

3.4 Polarization near the vortex

In the vicinity of the above-mentioned vortex, not only the surface brightness increases but the polarization degree changes as well. In particular, the polarization due to the emission of thermal radiation of aligned dust grains decreases with increasing surface brightness. This is shown in Fig. 8, which enlarges the section around the formed vortex. In contrast to the local environment, the dust density and thus the optical depth increase there. Consequently, the polarization due to the emission of thermal radiation and extinction of aligned dust grains cancel each other out. For the case of self-scattering, the resulting polarization vectors do not show a homogeneous pattern anymore due to the inhomogeneous radiation field. Moreover, the degree of polarization increases with increasing surface brightness.

Fig. 7 Wavelength dependence of the radial profile of the polarization degree. In the top and middle figure, only 100 μm and 1 mm-sized grains are aligned and contribute to the polarization, respectively. In the bottom figure, the polarization arises from self-scattered radiation only. The disk is inclined by 45∘. See Sects. 3.2 and 3.3 for details.

3.5 Flip of the polarization orientation

As shown in Fig. 4, the absorption efficiency along the minor grain axis becomes greater than the absorption efficiency along the major grain axis for particle sizes approximately the wavelength of observation. As a consequence, the resulting polarization orientation changes. In particular, for a wavelength of 3 mm and aligned grains with a size of 1 mm, the polarization orientation of the emitted radiation rotates by an angle of 90^∘ and is parallel to the minor axis of the disk, as shown in Fig. 9. Here, the polarization vectors coincide with the expected orientation of polarization due to self-scattering, as depicted in Figs. 6c and 6d. In contrast, if only 100 μm-sized grains are aligned, the orientation is still parallel to the major axis of the disk for all considered wavelengths. This is because the flip of the polarization is expected to occur only for wavelengths λ ≲ 2πa_eff.

Most interestingly, the polarization orientation also changes for self-scattered radiation at a wavelength of 3 mm (see Fig. 9). This polarization reversal (Yang et al. 2016) is also caused by dust grains, which are about the size of the observing wavelength. In particular, the polarization of scattered radiation is parallel to the scattering plane, causing rotation by an angle of 90^∘ in the polarization. However, the degree of linear polarization of thermally emitted radiation is similar or somewhat lower compared to the polarization due to self-scattering, that is, ≲ 1%. Therefore, if both mechanisms are considered at the same time, contributions of polarized emission and self-scattering may cancel each other out.

Fig. 8 Surface brightness map at a wavelength of 10 mm with an enlarged image section around the vortex. The polarization arises from thermal emission and dichroic absorption of aligned 1 mm-sized dust grains (top) and from self-scattered radiation (bottom) only. The disk is inclined by 45∘. See Sect. 3.4 for details.

Fig. 9 Surface brightness map of the total emitted thermal radiation at a wavelength of 3 mm overlaid with corresponding polarization vectors. The polarization arises from thermal emission and dichroic absorption of aligned 1 mm-sized dust grains (top) and from self-scattered radiation (bottom) only. The disk is inclined by 45∘. See Sect. 3.5 for details.

3.6 Comparison to an analytic drift model

Finally, we compare our results with those expected in the case of an analytic drift model. In contrast to numerical simulations, here the radial, azimuthal, and polar components of the headwind or drift velocity can be written as $$v_r=\frac{\mathrm{S}\mathrm{t}}{1+{\mathrm{S}\mathrm{t}}^2}\eta v_{\mathrm{K}},v_{\varphi }=-\frac{{\mathrm{S}\mathrm{t}}^2}{1+{\mathrm{S}\mathrm{t}}^2}\eta v_{\mathrm{K}},v_{\theta }\equiv 0\mathrm{m}{\mathrm{s}}^{-1},$$(9)

respectively, where St is the Stokes number and ηv_K is the gas rotation velocity relative to the Keplerian velocity (Kataoka et al. 2019). Here, we assume η v_K = 53 m s⁻¹. In this theoretical drift model, instabilities and turbulences such as the VSI are ignored.

However, the orientation of the drift velocity vector depends on the Stokes number, which causes the radial or azimuthal component to dominate (see Kataoka et al. 2019).

For the given total disk gas mass (see Table 1), Stokes numbers would be in the order of ≲ 0.1 for the largest 1 mm-sized dust grains. However, for these values, Eq. (5) is not satisfied. Thus, the grains are not aligned because of METs and the thermal radiation is unpolarized. Conversely, this means that with our disk model, a toroidal polarization pattern can only be produced if our initial assumptions for alignment (see Sect. 2.2) are not fulfilled. In this case, grain alignment due to RATs is an alternative explanation for the azimuthal orientation of polarization (Tazaki et al. 2017). In order to present exemplary results of polarization maps, we apply larger Stokes numbers of 1 and 10. For a Stokes number of 1, the polarization pattern is defined by a transition between a radial and azimuthal orientation. In contrast, a Stokes number of 10 causes a purely radial polarization pattern.

Figure 10 shows the surface brightness map of the emitted thermal radiation at a wavelength of 10 mm. The polarization is caused by elongated 1 mm-sized grains that are aligned to the drift velocity vector defined by Eq. (9).

The main difference from the results discussed in Sect. 3.2 is the orientation of the polarization. In particular, the analytic drift model assumes only radial and azimuthal drift velocity components, whereas the strength of each component depends on the Stokes number. This is in contrast to the results obtained from hydrodynamical simulations, since the VSI causes a dominating vertical drift velocity component that produces no radial or azimuthal polarization structure.

Moreover, the degree of linear polarization is significantly higher if the alignment is described by the analytic drift velocity model (Eq. (9)) compared to the polarization if the alignment is described by the velocity field obtained with the PLUTO code. This is due to the fact that the analytic velocity field described by Eq. (9) is homogeneous throughout the model space. In contrast, the velocity field from numerical hydrodynamical simulations is irregular with disturbances on small scales, causing a decrease in the net degree of polarization.

Fig. 10 Surface brightness map of the total emitted thermal radiation at a wavelength of 10 mm overlaid by polarization vectors assuming an analytic description for the drift velocity (see Eq. (9)). The Stokes number is set to St = 1 (top) and St = 10 (bottom). The inclination of the disk is 45∘. Here, only 1 mm-sized grains are aligned. See Sect. 3.6 for details.

3.7 Comparison to observations

Do we observe polarization due to aligned grains? With our 3D radiation hydrodynamical models, we obtain a detailed and rich dataset of the kinematic structure of gas and dust of a VSI turbulent protoplanetary disk. For a VSI turbulent disk, the drift velocity or headwind of the dust grains is dominated by the largescale vertical motions characteristic of the VSI. We found that in this model, if the grains are mechanically aligned, the polarization patters are mostly parallel to the major axis of the disk for a large range of millimeter wavelengths (see Figs. 6a and 6b).

However, in recent observations of the polarized emission of several protoplanetary disks at ALMA Band 7 (Stephens et al. 2017; Harrison et al. 2024; Lin et al. 2024), the polarization vectors are parallel to the projected minor axis of the protoplanetary disk, which is more similar to the self-scattering pattern (see Figs. 6c and 6d). At longer wavelengths, in particular at ALMA Band 3, the grains are aligned with their long axis along the azimuthal direction of the disk (Kataoka et al. 2017; Stephens et al. 2017), probably caused by RATs and not the mechanism considered in this work. Finally, for the disk around L1448 IRS3B, Looney et al. (2025) observed polarization that is consistent with the mechanical alignment of dust grains along the spiral dust structures.

More observations of further less massive disks and different instability simulations (Lesur et al. 2023) are needed to refine our understanding of how grains are aligned and how they emit. We also want to remind the reader that smaller grains that are less settled might shield the emission of larger grains (Sierra & Lizano 2020) and thus dominate the self-scatter emission. This could be the case of 100 μm-sized grains that dominate the self-scatter emission, shielding the emission of millimeter-sized grains drifting at the midplane.

4 Model assumptions and limitations

There are several limitations in our modeling of the mechanical alignment of large BAM dust aggregates that may affect the synthetic observations presented in this paper. First, the Monte Carlo nature of the MAD simulations introduces statistical noise in the angular velocity and alignment angle distributions³. We use 30 grain realizations per size bin, which we consider sufficient to capture the mean alignment behavior. This is because the Rayleigh reduction factor ℛ is computed using each of the 30 individual Monte Carlo realizations of Θ (see Eq. (4)). As a result, ℛ is an averaged quantity that smooths out both individual outliers in Θ and the expected Monte Carlo noise. Increasing the number of realizations would reduce the influence of this noise and could be explored in future studies that focus on the statistical properties of alignment.

Second, although our alignment simulations explicitly cover grain sizes up to 1.4 μm, we extrapolate the derived alignment parameters to grains up to 1 mm in size, assuming that alignment behavior converges for larger grains (see Appendix A). This assumption is motivated by the competing effects of increased surface area and a more spherical morphology, which tend to flatten the torque response across larger sizes. Although our approach follows the trends reported in Reissl et al. (2023), this extrapolation remains an assumption and should be tested directly in future work through high-resolution MAD simulations of larger aggregates.

Finally, our radiative transfer modeling is based on a single hydrodynamical snapshot taken after 750 orbits of disk evolution. This choice is justified by the assumption that the alignment timescale is shorter than the local dynamical timescales. However, a more detailed comparison between the alignment timescale and the temporal variability of the disk environment would strengthen this justification. In principle, time-averaged snapshots may offer a more robust basis for post-processing if temporal variability proves significant. Although our approach is consistent with previous studies, a more comprehensive, timeresolved treatment of grain alignment may be desirable in future investigations.

We emphasize once again that our modeling of dust polarization focuses exclusively on two distinct effects: (i) dust scattering and (ii) polarized thermal emission from mechanically aligned grains. However, it is highly unlikely that dust polarization in protoplanetary disks arises from one single mechanism in isolation. Multiple alignment and polarization mechanisms have been proposed over the years and are expected to operate simultaneously within disk environments (see, e.g., Andersson 2015). Consequently, a comprehensive model of the dust polarization process would require incorporating all relevant mechanisms into a unified RT framework (Hoang & Lazarian 2016; Giang et al. 2023). However, the question of how all these mechanisms work in tandem is still a matter of debate. Hence, such an endeavor is beyond the scope of the present study and will be addressed in future work.

5 Summary and conclusions

The goal of this study is to evaluate the impact of mechanical grain alignment in protoplanetary disks on the potentially observable polarization. The grain alignment is assumed to be the result of the drift motion of the grains. For this purpose, we combined radiation hydrodynamical simulations to determine the gas and dust density distribution as well as their velocity fields, Monte Carlo dust-gas interaction simulations to calculate the mechanical alignment of dust in a gas flow, and finally Monte Carlo polarized radiative transfer simulations to obtain synthetic observations for the flux and its degree of polarization of thermally re-emitted radiation of aligned nonspherical grains as well as radiation scattered by dust grains.

First, we applied the code PLUTO, which solves the radiation hydrodynamic equations and provides the density distribution (Fig. 3) and the velocity components (see Fig. C.1). Based on these quantities, the headwind or drift velocity was calculated (see Eq. (1)), which defines the direction of grain alignment due to mechanical alignment. Subsequently, by applying ballistic aggregation and migration of the primary particle, dust grain aggregates were created (see Sect. 2.1). Next, Monte Carlo dust-gas interaction simulations were used to determine a lower and upper grain size limit for alignment in a gas flow (see Sect. 2.2). Moreover, we estimated the effect of depolarization due to grain precession via the Rayleigh reduction factor (see Eq. (5)). Finally, POLARIS calculated the temperature of dust grains in the protoplanetary disk and then the thermal and self-scattered radiation by applying Monte Carlo radiative transfer simulations (see Sect. 2.4). For the calculation of polarization due to dichroic emission and absorption, nonspherical grains were represented by oblate shaped particles that are aligned with their minor axis to the drift velocity vector. The optical properties of these dust grains were calculated using DDSCAT, whereas for temperature and scattering simulations, the optical properties were calculated using the Mie scattering theory (see Sect. 2.3).

We find that large dust grains, which contribute the most to the observable polarization, fulfill the alignment criteria (see Eq. (5) and Fig. 5) and are therefore aligned due to the gas flow. In particular, because of hydrodynamical instabilities such as the VSI, the polar drift velocity that is perpendicular to the disk midplane is the dominant velocity component (see Fig. C.1). As a result, most grains align with their short axis parallel to the rotational axis of the disk. Consequently, the orientation of the polarization of the thermal emission is close to zero for a face-on disk, but parallel to the major disk axis for an inclined disk (see Fig. 6). Here, the highest degree of polarization is found for the longest considered wavelength of 10 mm. Most notably, the orientation of polarization is in contrast to analytic drift velocity models, which assume solely radial and azimuthal drift components (see Sect. 3.6).

In addition, polarization arises from self-scattering as well (see Sect. 3.3). In this case, the orientation of the polarization is usually parallel to the projected minor disk axis. Here, the contribution of self-scattering is largest for the shortest considered wavelength of 0.87 mm. With increasing wavelength, and therefore decreasing optical depth, the impact of self-scattering becomes less important in the outer disk regions.

Moreover, the VSI potentially forms vortices that can trap particles. Due to the increase in dust density, the surface brightness increases near the vortex as well (see Sect. 3.4). In contrast, the degree of polarization from thermally emitted radiation decreases, whereas the polarization from self-scattering produces a nonuniform pattern because of the anisotropic radiation field.

Finally, we observe a change in the orientation of the polarization by an angle of 90^∘. This behavior is expected if the observing wavelength is about the size of the dust grains, that is, λ ≲ 2 π a. In particular, at a wavelength of 3 mm, the orientation changes in the polarization of thermal emission of aligned 1 mm-sized dust grains as well as in the polarization due to self-scattering (see Fig. 9). Most interestingly, the polarization vectors for the thermally emitted radiation are parallel to the projected minor disk axis and coincide with the expected case of polarization due to self-scattering, as depicted in Fig. 6.

In conclusion, grain alignment driven solely by mechanical torques proves to be an efficient mechanism for aligning millimeter-sized grains in protoplanetary disks, particularly when hydrodynamical instabilities such as the VSI are taken into account. Our synthetic polarimetric observations in the millimeter wavelength range demonstrate that the resulting polarization may serve as a promising diagnostic tool to trace the orientation of aligned millimeter-sized dust grains.

Acknowledgements

The authors thank the anonymous referee for useful and constructive comments. M.L.-S. and S.W. acknowledge financial support from the DLR/BMWK grant 50 OR 2210. S.R. acknowledge financial support from the Heidelberg cluster of excellence (EXC 2181 – 390900948) “STRUCTURES: a unifying approach to emergent phenomena in the physical world, mathematics, and complex data”, specifically via the exploratory project EP 4.4. Furthermore, S.R. thanks for funding from the European Research Council in the ERC synergy grant “ECOGAL – Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets” (project ID 855130). This research made use of Astropy (Astropy Collaboration 2013, 2018, 2022), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), scikit-image (van der Walt et al. 2014), and the Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC21M00561.

Appendix A Characteristic alignment timescale

In this section, we briefly evaluate the extrapolation of the alignment process for large grains, that is, aggregates with an effective radius of 100 μm and 1 mm, respectively.

The time-dependent magnitude of the angular velocity, ω(t), of each grain can be evaluated via $$\omega (t)={\omega }_{max}[1-\exp (-\frac{t}{{\tau }_{\text{drag\hspace{0.17em}}}})]$$(A.1)

by solving the differential equation Eq. (3). For simplicity, we use the maximum mechanical torque Γ_MET(Θ) over all angles Θ to estimate the upper limit of the maximum angular momentum, ω_max = max_Θ(Γ_MET(Θ)) τ_drag/I_m, of each individual BAM1 dust aggregate. Consequently, each aggregate reaches its maximum angular momentum, ω_max, within a 99% limit over a characteristic drag timescale of 5 × τ_drag.

In Fig. A.1, we present the ratio of the maximum angular velocity, ω_max, to the thermal angular velocity ω_th, averaged over all ensembles of dust aggregates with the same effective radius a_eff, for velocity ratios of v_drift/v_th = 0.1 and v_drift/v_th = 1.0. Although some ensembles may contain outliers, the ensemble average shows a clear trend in which the ratio ω_max/ω_th approaches an asymptotic value with increasing a_eff. This asymptotic behavior is consistent with the findings of Reissl et al. (2023), where the mechanical spin-up process was modeled using fractal dust aggregates. However, in Reissl et al. (2023), only aggregates with an upper grain radius of a_eff = 800 nm were considered, whereas in this publication we use aggregates up to a_eff = 1.4 μm. Using aggregates up to 1.4 μm reveals the asymptotic trend even more clearly.

We again emphasize that the simulation of mechanical alignment parameters beyond 1.4 μm aggregates is numerically unfeasible (see Reissl et al. 2023). To quantify the trend, we perform a least-squares fit with an arctan function. As shown in Fig. A.1, the ratio ω_max/ω_th already approximates its asymptote for aggregates with a_eff ≥ 1.0 μm. Hence, in this paper, we assume that grains with a_eff = 100 μm and 1 mm have a similar ratio of ω_max/ω_th as the much smaller aggregates within the range a_eff = 1.0 μm to 1.4 μm.

We extrapolate the characteristic drag timescale, τ_drag, for the considered aggregate sizes of 100 μm and 1 mm similarly. Since the drag depends only marginally on IR photon emission due to the low dust temperature within the disk, it is dominated by the gas drag timescale ${\tau }_{\text{gas\hspace{0.17em}}}\propto I_{\mathrm{m}}{n}_{\mathrm{g}}^{-1}{a}_{\text{eff\hspace{0.17em}}}^{-4}$ (see e.g. Draine & Weingartner 1996; Jäger et al. 2024). Given that the maximal moment of inertia of each aggregate scales with $I_{\mathrm{m}}\propto a_{\text{eff\hspace{0.17em}}}^5$, we obtain a nearly linear dependence of τ_drag ∝ a_eff.

In Fig. A.1, we show τ_drag as a function of a_eff for the ratios v_drift/v_th = 0.1 and v_drift/v_th = 1.0, respectively. Based on the average τ_drag across all aggregate ensembles of similar size, we linearly extrapolate τ_drag up to the utilized grain radii of a_eff = 100 μm and 1 mm, respectively. For a ratio of v_drift/v_th = 0.1, the large aggregates exhibit drag timescales of 1.71 × 10⁶ s and 1.70 × 10⁷ s⁴. In contrast, for v_drift/v_th = 1.0, the drag timescales decrease to 1.54 × 10⁵ s and 1.53 × 10⁶ s, respectively, assuming a gas number density of n_g = 10¹⁰ cm⁻³ and a gas temperature of T_g = 20 K. Hence, the largest considered grain size with a_eff = 1 mm would reach its final grain alignment configuration, that is, an angular velocity ω_max, within 5 × τ_drag ≈ 3 yr at most based on the extrapolated data from MAD simulations. For increasing gas density, these timescales become even shorter. We emphasize that micro-turbulence of the gas on aggregate scales is already accounted for in the MAD model, where the magnitude and direction of the gas velocity are sampled from a distribution function based on gas temperature and drift velocity.

Therefore, even the largest aggregates reach their final alignment configuration on timescales significantly shorter than that of the vertical shear instability, which operates on the order of tens of orbital timescales, which represents several thousands of years at 50 au.

Fig. A.1 Average alignment over all BAM1 grain shapes according to the MAD model (dot shapes) and their extrapolation (diamond shapes) for the drag timescale tdrag (red) and the ratio of angular velocity ω (blue) with respect to the thermal angular velocity ωth dependent on the effective grain radius aeff. The data is modeled for the fiducial values of ng = 1010 cm−3, Tg = 20 K, and Td = 20 K. The left panel depicts data for the distinct ratio of the drift velocity vdrift to the thermal velocity vth with vdrift/vth = 0.1 while the right panel if for vdrift/vth = 1.0. The best fits (dashed lines) of τdrag are performed via τdrag = a ×(aeff/1 nm)+b and ω/ωth = a × arctan [b × log10(aeff/1 nm)+c]+d, respectively. The MAD data is extrapolated data for the two particular grain species of aeff = 100 μm as well as aeff = 1 mm utilized in this paper.

Appendix B Direction of grain alignment

Usually, grain alignment and the resulting dust polarization are estimated using toy models or dust grains with one idealized shape throughout (Purcell & Spitzer 1971; Das & Weingartner 2016; Hoang et al. 2018). The advantage of the MAD framework is that it considers the mechanical alignment of a large ensemble of irregularly shaped grains, enabling the evaluation of the expected dust polarization based on the resulting distribution of the alignment angle Θ.

In Fig. B.1 we present the distribution of Θ for selected grain ensembles with a distinct effective radius between a_eff = 100 nm and 1.4 μm dependent on the two different velocity ratios of v_drift/v_th = 0.1 and v_drift/v_th = 1.0. For the ratio v_drift/v_th = 0.1 and grains with a_eff = 100 nm, the distribution is very narrow close to Θ=0^∘. Hence, the resulting Rayleigh reduction factor amounts to ℛ=0.692. For a ratio of v_drift/v_th = 1.0, the Rayleigh reduction factor decreases to ℛ=0.311. In this particular case, the lower value of ℛ is because part of the aggregates within the ensemble becomes rotationally disrupted (compare Fig. 2). Larger aggregates with a_eff ≥ 1 μm rotate slower in our model because they are more spherically shaped and cannot be rotationally destroyed for v_drift/v_th ≲ 10. However, as shown in Fig. B.1 the distribution of Θ for aggregates with a_eff ≥ 1 μm becomes even broader compared to smaller grains with a second peak appearing closer to Θ=90^∘. Here, the resulting Rayleigh reduction factor is even lower with an average value of ℛ= 0.29. We report similar distributions of Θ for the entire range 1.78 × 10⁻² ≲ v_drift/v_th ≲ 10.

Although we observe convergence in the Rayleigh reduction factor, ℛ, for grain sizes larger than approximately 1 μm, we note that individual angle distributions can exhibit bimodal features, particularly under near-sonic drift conditions, as shown in Fig. B.1. These features resemble the classical Gold-type alignment, in which grains align with their short axis perpendicular to the drift velocity (Gold 1952a,b). However, our results show a much broader distribution across the full range 0^∘ < Θ < 90^∘, including a peak near Θ=0^∘, which is not captured by the classical Gold mechanism. In this regard, the MAD framework provides an improvement by accounting for more realistic, irregular grain shapes and their complex gas-grain interactions, even compared to more recent mechanical alignment models (see, e.g., Lazarian & Hoang 2007b; Das & Weingartner 2016; Hoang et al. 2018). Moreover, the drift velocities in the utilized MHD snapshot are around 0.1 to 0.2 Mach, that is, subsonic i.e. Gold alignment would not predict any dust polarization at all.

In summary, we see a convergence of angular momentum for larger grains (see Figs. 2 and A.1) and characteristic grain alignment timescales that are significantly shorter than those of the VSI. Consequently, for the utilized grain sizes of 100 μm and 1 mm, we extrapolate a Rayleigh reduction factor of ℛ ≈ 0.29 within the velocity ratio range 1.78 × 10⁻² ≲ v_drift/v_th ≲ 10.

Fig. B.1 Distribution of the alignment angle Θ (compare Fig. 1) of a distinct ensemble of BAM1 dust aggregates according to the MAD framework. The four rows (top, second, third, bottom) show the distribution for an effective radius of aeff = 100 nm, aeff = 1000 nm, aeff = 1200 nm, and aeff = 1400 nm, respectively. Individual alignment angles Θ are simulated for fiducial values of ng = 1010 cm−3, Tg = 20 K, and Td = 20 K whereas the left column is for the ratio of velocities vdrift/vth = 0.1 while the right column is for the ratio vdrift/vth = 1.0. The Rayleigh reduction factor ℛ (see Eq. 4) is calculated considering both the distribution of Θ within 0∘ to 90∘ and the ratio ζ of grains that either do not align with the velocity field at all because of slow rotation or grains that are rotationally destroyed.

Appendix C Velocity field

Figure C.1 shows the three components of the drift velocity vector in the disk midplane. Here, the polar drift velocity is the dominant component, causing alignment with the short axis of the dust grains perpendicular to the disk midplane. The alternating polar velocity of the 1 mm-sized grains at the inner disk rim is a numerical effect caused by the boundary conditions.

Fig. C.1 Radial (left), azimuthal (middle), and polar (right) drift velocity component in the disk midplane of 100 μm (top) and 1 mm-sized (bottom) grains. Note the different scales of the color bar for each component of the drift velocity. The gray dashed line indicates the radius of the disk at 100 au.

Appendix D Dust temperature distribution

Figure D.1 shows the average dust grain temperature of the circumstellar disk model in a cross-section view resulting from a simulation using POLARIS. The temperature values are averaged azimuthally. At the inner rim of the disk, the temperatures are up to about 100 K.

Fig. D.1 Cross-section view of the average dust grain temperature. The temperature values are averaged azimuthally.

References

All Tables

All Figures

Fig. 1 Exemplary BAM1 dust aggregate with an effective radius of aeff = 1.2 μm. The aggregate rotates with an angular velocity, ω, that corresponds to its maximum moment of inertia, Im. The rotation and precession of the dust aggregate result from gas-dust interactions, where the gas and dust components move with a relative velocity, vdrift, through the disk. The angle, Θ, is defined as the angle between ω and vdrift.

In the text

Fig. 2 Left panel: Angular velocity, ω, as a function of the ratio vdrift/vth. The effective grain radii, aeff = 0.1 μm to 1.4 μ m, are color-coded. Solid lines represent the average values of ω, while the shaded areas indicate the range between the minimum and maximum values. Dotted lines show the lower limit, 3 × ωth, for long-term stable alignment, and dashed lines represent the upper limit for rotational disruption, ωdisr. Right panel: Same as the left panel, but for the Rayleigh reduction factor ℛ (solid lines). Dashed-dotted lines depict the parametric representation of ℛ for grains with aeff ≥ 1 μm. See Sect. 2.2 for details.

In the text

	Fig. 3 Cross-section view of the dust mass density in units of g cm−3 of 0.1 μm to 10 μ m-sized grains (green), 100 μ m-sized grains (purple), and 1 mm-sized grains (orange). The values are azimuthally averaged.
In the text

Fig. 4 Ratio of absorption efficiencies Q⊥/Q‖(top) and absorption efficiency (bottom) of an oblate spheroid as a function of radius of an equal volume sphere. This figure shows an axis ratio of the spheroid of 1/1.5 and a wavelength of radiation of 3 mm. The absorption efficiency Q‖is the efficiency along the minor grain axis, that is, the axis of alignment. If a/λ ≥ 1, we assume a spherical particle, thus, Q⊥ = Q‖. See Sect. 2.3 for details.

In the text

	Fig. 5 Ratio of the drift velocity and the thermal velocity in the disk midplane of 100 μm (top) and 1 mm-sized (bottom) grains. Dust grains are aligned if 1.78 × 10−2 ≤ vdrift/vth ≤ 10 (see Eq. (5)), otherwise they are marked in red in the figure. See Sect. 3.1 for details.
In the text

Fig. 6 Surface brightness map of the total emitted thermal radiation at a wavelength of 1.3 mm (left) and 10 mm (right) overlaid with polarization vectors at these wavelengths. In the top row (a, b), the polarization arises from thermal emission and dichroic extinction of aligned dust grains only. At wavelengths of 1.3 mm and 10 mm, only 100 μm and 1 mm-sized grains are assumed to be aligned, respectively. In contrast, in the bottom row (c, d), the polarization arises from self-scattered radiation only. The disk is inclined by 45∘. See Sects. 3.2 and 3.3 for details.

In the text

	Fig. 7 Wavelength dependence of the radial profile of the polarization degree. In the top and middle figure, only 100 μm and 1 mm-sized grains are aligned and contribute to the polarization, respectively. In the bottom figure, the polarization arises from self-scattered radiation only. The disk is inclined by 45∘. See Sects. 3.2 and 3.3 for details.
In the text

	Fig. 8 Surface brightness map at a wavelength of 10 mm with an enlarged image section around the vortex. The polarization arises from thermal emission and dichroic absorption of aligned 1 mm-sized dust grains (top) and from self-scattered radiation (bottom) only. The disk is inclined by 45∘. See Sect. 3.4 for details.
In the text

	Fig. 9 Surface brightness map of the total emitted thermal radiation at a wavelength of 3 mm overlaid with corresponding polarization vectors. The polarization arises from thermal emission and dichroic absorption of aligned 1 mm-sized dust grains (top) and from self-scattered radiation (bottom) only. The disk is inclined by 45∘. See Sect. 3.5 for details.
In the text

	Fig. 10 Surface brightness map of the total emitted thermal radiation at a wavelength of 10 mm overlaid by polarization vectors assuming an analytic description for the drift velocity (see Eq. (9)). The Stokes number is set to St = 1 (top) and St = 10 (bottom). The inclination of the disk is 45∘. Here, only 1 mm-sized grains are aligned. See Sect. 3.6 for details.
In the text

Fig. A.1 Average alignment over all BAM1 grain shapes according to the MAD model (dot shapes) and their extrapolation (diamond shapes) for the drag timescale tdrag (red) and the ratio of angular velocity ω (blue) with respect to the thermal angular velocity ωth dependent on the effective grain radius aeff. The data is modeled for the fiducial values of ng = 1010 cm−3, Tg = 20 K, and Td = 20 K. The left panel depicts data for the distinct ratio of the drift velocity vdrift to the thermal velocity vth with vdrift/vth = 0.1 while the right panel if for vdrift/vth = 1.0. The best fits (dashed lines) of τdrag are performed via τdrag = a ×(aeff/1 nm)+b and ω/ωth = a × arctan [b × log10(aeff/1 nm)+c]+d, respectively. The MAD data is extrapolated data for the two particular grain species of aeff = 100 μm as well as aeff = 1 mm utilized in this paper.

In the text

Fig. B.1 Distribution of the alignment angle Θ (compare Fig. 1) of a distinct ensemble of BAM1 dust aggregates according to the MAD framework. The four rows (top, second, third, bottom) show the distribution for an effective radius of aeff = 100 nm, aeff = 1000 nm, aeff = 1200 nm, and aeff = 1400 nm, respectively. Individual alignment angles Θ are simulated for fiducial values of ng = 1010 cm−3, Tg = 20 K, and Td = 20 K whereas the left column is for the ratio of velocities vdrift/vth = 0.1 while the right column is for the ratio vdrift/vth = 1.0. The Rayleigh reduction factor ℛ (see Eq. 4) is calculated considering both the distribution of Θ within 0∘ to 90∘ and the ratio ζ of grains that either do not align with the velocity field at all because of slow rotation or grains that are rotationally destroyed.

In the text

	Fig. C.1 Radial (left), azimuthal (middle), and polar (right) drift velocity component in the disk midplane of 100 μm (top) and 1 mm-sized (bottom) grains. Note the different scales of the color bar for each component of the drift velocity. The gray dashed line indicates the radius of the disk at 100 au.
In the text

	Fig. D.1 Cross-section view of the average dust grain temperature. The temperature values are averaged azimuthally.
In the text