© The Authors 2025

1. Introduction

Many astrophysical systems consist of a central magnetized star surrounded by an accretion disk. One of the most widely accepted magnetospheric accretion models to date is that of Ghosh & Lamb (1978), which assumes that near the star, magnetic fields are so strong that the plasma is forced to corotate with the star. The model relies on the balance between the spin-up torque on the star, which is due to the interaction with the region inside the corotation radius, and the spin-down torque from the region outside it. It includes field shearing due to the disk flow. In this model, the growth of the toroidal field component (B_ϕ) is limited by reconnection, and it is assumed that a steady state is reached. Many of the predicted features have been confirmed through numerical simulations, and new insights into the accretion process onto a star with a dipole field have been revealed. For example, it was found that matter can fall directly onto the star’s equator (Mishra et al. 2023) or form an accretion column at a higher latitude, the so-called funnel effect (Romanova et al. 2002; Zanni & Ferreira 2009; Čemeljić 2019), but it can also escape from the vicinity of the star in the form of a stellar wind or a conical and axial outflow (Zanni & Ferreira 2013; Kotek et al. 2020; Čemeljić & Brun 2023).

However, several observational results indicate that the configuration of the magnetic field of strongly magnetized stars may not be dipolar. For instance, Bouvier et al. (2006) argued that there are many signs of magnetospheric accretion in classical T Tauri stars and that the field is probably dipolar at larger distances from the star, but may have a strong multipolar component close to the star. Donati et al. (2007) derived the magnetic topology on the surface of the classical T Tauri star V2129 Oph, and found that it has a dominant octupole component with B_oc ≈ 1.2 kG and a weaker dipole component with B_d ≈ 0.35 kG. Additionally, early research on X-ray emission in neutron stars also suggested a quadrupole magnetic component (Shakura et al. 1991; Panchenko & Postnov 1994). Thus, the magnetic field configuration likely differs from star to star depending on the significance of the dipole component. Lipunov (1978) suggested a framework to extend the spherical accretion model to the case of stellar magnetic fields with higher-order multipole moments. von Rekowski & Brandenburg (2006) performed axisymmetric simulations considering the disk–magnetosphere interaction when the magnetic field is generated by the dynamo effect and found a time-varying stellar magnetic field with a complex multipolar configuration.

Although observational and theoretical evidence suggests that stellar magnetic configurations may not be purely dipolar, most simulations of disk–magnetosphere interactions assume a dipolar field configuration. Exceptions include the nonrelativistic simulations of Long et al. (2007, 2008) and Cieciuch & Čemeljić (2022) and the general relativistic simulations of Das et al. (2022). Çıkıntoğlu (2023) studied the interaction between the disk and a neutron star with a quadrupole field configuration, concluding that the effect of the quadrupole depends on the location of the inner radius of the disk and the strength of the quadrupole field.

In this work, we examined the magnetospheric accretion through 2.5-D magnetohydrodynamic (MHD) simulations, exploring various initial magnetic configurations. These include dipolar (D), quadrupolar (Q), octupolar (O) fields, as well as combined dipolar-quadrupolar (DQ) and dipolar-octupolar (DO) field configurations, all within the context of disk accretion onto a rotating star. We considered the observational parameters for the star IRAS-17449+2320, which is characterized by a strong magnetic field (≈6.2 kG) and a low rotational velocity (Korčáková et al. 2022). This particular object has raised the possibility of post-mergers occurring among FS CMa stars, a subgroup of B[e] stars whose distinctive spectral characteristic is the existence of forbidden emission lines and a strong IR excess.

In our earlier work, we performed 2.5-D MHD simulations of disk accretion onto stars with pure dipole fields (Moranchel-Basurto et al. 2023, 2024), revealing distinct features in the funnel streams and associated gaps compared to cases involving a pure dipole field and magnetic diffusivity. In scenarios with significant resistivity, accretion always occurs along the midplane and produces magnetospheric ejections from the disk and stellar surface.

The paper is laid out as follows. In Section 2, we describe our numerical model and the different magnetic configurations considered. The results of the simulations are given in Section 3. In Section 4, we give a brief discussion. Finally, the conclusions are given in Section 5.

2. Model setup and magnetic field configurations

We modeled the evolution of the gas in a thin accretion disk and the magnetosphere region surrounding an FS CMa type star. The simulations solve the MHD equations (Eq. (1) to Eq. (4)), including both resistivity and viscosity,

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathit{\nabla }·(\rho \mathit{v})=0,\end{array}$$(1)

$$\begin{array}{c}\hfill \frac{\partial \rho \mathit{v}}{\partial t}+\mathit{\nabla }·(\rho \mathit{v}{\mathit{v}}^T-\frac{\mathit{B}{\mathit{B}}^T}{4\pi })+\mathit{\nabla }(P+\frac{B^2}{8\pi })+\rho \mathit{\nabla }\mathrm{\Phi }=0,\end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial e}{\partial t}& +\mathit{\nabla }·[\mathit{v}(e+P+\frac{B^2}{8\pi })-\frac{1}{4\pi }(\mathit{v}·\mathit{B})\mathit{B}+(\eta ·\mathit{J})\times \mathit{B}]\hfill \\ \hfill & =-\rho (\mathit{\nabla }\mathrm{\Phi })·\mathit{v},\hfill \end{array}$$(3)

and the induction equation

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}-\mathit{\nabla }\times (\mathit{v}\times \mathit{B}-\eta ·\mathit{J})=0,\end{array}$$(4)

where P is the thermal gas pressure, ρ the gas density, v the gas velocity vector, B the magnetic flux density vector, J the electric current density vector, η the electrical resistivity, which is related with the magnetic diffusivity defined as ν = η/4π, Φ = −GM_⋆/R the gravitational potential, and e the total energy density given by e = P/(γ − 1)+ρv²/2 + B²/(8π), where γ = 5/3 is the plasma polytropic index. The viscous stress tensor τ is defined as:

$$\begin{array}{c}\hfill \tau =\eta [(\mathrm{\nabla }v)+{(\mathrm{\nabla }v)}^T-\frac{2}{3}(\mathrm{\nabla }v)I]\end{array}$$(5)

where I is the unit tensor.

We employed a two-dimensional computational domain (R, θ) with three-dimensional vector fields in a spherical coordinate system (R, θ, ϕ), assuming axisymmetry around the stellar rotation axis, the so called 2.5-D approach. Our calculations were carried out with the code PLUTO¹ (Mignone et al. 2007), using the HLL Riemann solver and the second-order Runge–Kutta integration in time to advance conserved variables. To compute viscous and resistive terms, a second-order finite difference approximation, integrated explicitly in time, was adopted.

In Table 1 we present the values of parameters used in this work, and is divided into two sections. The first section presents the physical parameters adopted in the simulations for each component involved (the central star, the disk, and the corona). Since our study focuses on the analysis of the surroundings of FS CMa stars, the values considered for the central star are mainly based on the observational parameters for the object IRAS-17449+2320 reported by Korčáková et al. (2022). The second section of Table 1 provides a summary of the computational grid and the used normalization values, as the simulations are performed using dimensionless equations².

2.1. Model description

The density (ρ_d) and pressure (P_d) of the gas in the accretion disk are given by the following expressions (e.g., Zanni & Ferreira 2009; Moranchel-Basurto et al. 2024):

$$\begin{array}{c}\hfill {\rho }_d(R,r)={\rho }_{d0}{\left\{\frac{2}{5h^2}[\frac{R_0}{R}-(1-\frac{5h^2}{2})\frac{R_0}{r}]\right\}}^{3/2},\end{array}$$(6)

$$\begin{array}{c}\hfill P_d=h^2{\rho }_{d0}{v}_{K0}^2{\left(\frac{{\rho }_d}{{\rho }_{d0}}\right)}^{5/3},\end{array}$$(7)

where ρ_d0 and v_K0 are the density and Keplerian velocity at the midplane of the disk at R = R₀, h = C_s/v_K is the disk aspect ratio, with C_s the isothermal sound speed and $v_K=\sqrt{G{M}_{\star }/r}$. The values of these parameters are given in Table 1. Note that r = R sin θ is the cylindrical radius.

The initial disk atmosphere is defined with the assumption of a polytropic gas in hydrostatic equilibrium, without rotation. The density and pressure are given by,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{atm}}(R)={\rho }_{\mathrm{atm}}^0{\left(\frac{R_{\star }}{R}\right)}^{\frac{1}{\gamma -1}},\end{array}$$(8)

and

$$\begin{array}{c}\hfill P_{\mathrm{atm}}(R)={\rho }_{\mathrm{atm}}^0\frac{\gamma -1}{\gamma }\frac{G{M}_{\star }}{R_{\star }}{\left(\frac{R_{\star }}{R}\right)}^{\frac{\gamma }{\gamma -1}},\end{array}$$(9)

with γ = 5/3 (e.g., Zanni & Ferreira 2009). The ratio of the coronal to disk density at R = R₀ and z = 0 is ρ_atm⁰/ρ_d0 = 0.01, a value adopted in all models.

To model the star’s rotation and the boundary conditions of a perfect accretor, we followed the descriptions given in previous works (Zanni & Ferreira 2009; Čemeljić 2019; Moranchel-Basurto et al. 2024).

2.2. Magnetic configurations

In our work, the magnetic configurations exhibit axial symmetry around the rotational axis of the star and disk. Additionally, we assumed that the magnetic moments are aligned with the stellar rotational axis, which allowed us to express the magnetic field components in terms of the polar magnetic field strength (see Gregory et al. 2010 for a derivation) as follows.

For the pure dipolar case:

$$\begin{array}{c}\hfill B_R^d(R,\theta )=B_{\star }^d{\left(\frac{R_{\star }}{R}\right)}^{3}cos\theta ,\end{array}$$(10)

and

$$\begin{array}{c}\hfill B_{\theta }^d=\frac{1}{2}{B}_{\star }^d{\left(\frac{R_{\star }}{R}\right)}^{3}sin\theta .\end{array}$$(11)

For a quadrupolar configuration, the components are:

$$\begin{array}{c}\hfill B_R^q(R,\theta )=\frac{1}{2}{B}_{\star }^q{\left(\frac{R_{\star }}{R}\right)}^4(3{cos}^2\theta -1),\end{array}$$(12)

and

$$\begin{array}{c}\hfill B_{\theta }^q=B_{\star }^q{\left(\frac{R_{\star }}{R}\right)}^{4}sin\theta cos\theta .\end{array}$$(13)

Finally, for an octupolar field configuration:

$$\begin{array}{c}\hfill B_R^{\mathit{oc}}(R,\theta )=\frac{1}{2}{B}_{\star }^{\mathit{oc}}{\left(\frac{R_{\star }}{R}\right)}^5(5{cos}^3\theta -3cos\theta ),\end{array}$$(14)

and

$$\begin{array}{c}\hfill B_{\theta }^{\mathit{oc}}(R,\theta )=\frac{3}{8}{B}_{\star }^{\mathit{oc}}{\left(\frac{R_{\star }}{R}\right)}^5(5{cos}^2\theta sin\theta -sin\theta ).\end{array}$$(15)

In Equations (10)–(15), B_⋆^d, B_⋆^q, and B_⋆^oc are the polar magnetic field strengths of the dipole, quadrupole, and octupole at the stellar rotation pole, respectively.

Additionally, we explored models that combine a dipole with a quadrupole configuration, as well as a dipole with an octupole configuration. More specifically, the magnetic field is given by B^dq = B^d + B^q for dipole plus quadrupole, and B^do = B^d + B^oc, for the dipole plus octupole configuration.

A compilation of the potential vectors (A) and the magnetic flux functions (ψ) for the dipolar, quadrupolar and octupolar configurations is given in Appendix A.

To specify the total magnetic field strength on the stellar surface, we used a dimensionless parameter μ, defined as the ratio of the measured magnetic field strength on the stellar surface (B_⋆) to the reference field strength $B_0\equiv \sqrt{4\pi {\rho }_{d0}{V}_{K0}^2}$. For B_⋆ = 6.2 kG and for the values of ρ_d0 and V_K0 adopted in this work (see Table 1), B₀ = 69.42 G and therefore μ = B_⋆/B₀ = 89.3.

2.3. Diagnostics

2.3.1. Angular momentum flux

To analyze the different behaviors in the mass flow accreted by the star through this study, we considered the angular momentum fluxes carried by the magnetic field f_B, and by the matter f_m, as given by

$$\begin{array}{c}\hfill {\mathbf{f}}_B=\frac{r{B}_{\varphi }{\mathbf{B}}_p}{4\pi }\end{array}$$(16)

and

$$\begin{array}{c}\hfill {\mathbf{f}}_m=-\rho r{v}_{\varphi }{\mathbf{v}}_p,\end{array}$$(17)

respectively, with B_p and v_p the poloidal magnetic field and poloidal velocity vectors.

2.3.2. Disk truncation radius

The most commonly adopted definition of the truncation radius R_trunc is inferred from the condition: $B^2/(8\pi )={\rho }_d{\mathit{v}}_K^2$, with B and ρ_d evaluated at the midplane of the disk (e.g., Hartmann et al. (2016), Zhu et al. (2024)). If the magnetic field is a dipole,

$$\begin{array}{c}\hfill \frac{R_{\mathrm{trunc},\mathrm{th}}}{R_{\star }}={\left(\frac{{(B_{\star }^d)}^2{R}_{\star }}{32\pi G{M}_{\star }{\rho }_{d0}}\right)}^{2/7},\end{array}$$(18)

whereas for a pure octupole magnetic field,

$$\begin{array}{c}\hfill \frac{R_{\mathrm{trunc},\mathrm{th}}}{R_{\star }}={\left(\frac{9{(B_{\star }^{\mathit{oc}})}^2{R}_{\star }}{2^8\pi G{M}_{\star }{\rho }_{d0}}\right)}^{2/15}.\end{array}$$(19)

Using the values provided in Table 1, if the magnetic field is dipolar with a strength of $B_{\star }^d=6.2$ kG, the truncation radius is R_trunc = 3.48R_⋆. Assuming instead an octupolar magnetic field of the same strength, we obtain R_trunc = 1.82R_⋆.

We emphasize that the initial inner radius of the disk, R_d, for the models with configurations dipolar, quadrupolar and the compositions of dipolar plus quadrupolar field were fixed initially at R = 3.48R_⋆, which is the value obtained for the truncation radius considering a purely dipolar magnetic configuration using Eq. (18) and the values in Table 1. It is also equivalent to ∼0.75 times the corotation radius, $R_{\mathrm{co}}={(G{M}_{\star }/{\mathrm{\Omega }}_{\star }^2)}^{\frac{1}{3}}=4.65$, calculated from the same Table 1. On the other hand, for the models with dipolar plus octupolar configuration we considered that the inner radius of the disk is located at R = 1.82R_⋆.

3. Results

In this section we present the results of the different models labeled in Table 2. Figure 1 shows the initial gas density (in units of ρ₀) and the initial magnetic flux lines (white lines) for five different magnetic strength ratios. This figure shows each of the different shapes followed by the magnetic field lines, based on the contribution of each magnetic moment.

Fig. 1. Initial density (in units of ρ0) and magnetic field lines (white curves) for model D (first panel), model Q (second), model O (third), model DQ (fourth) and model DO (fifth). In all cases, the sum of the different multipolar contributions to the stellar magnetic field is equal to 89.3. Here, r and z represent the cylindrical coordinates, measured in units of the stellar radius.

3.1. Pure dipole and pure quadrupole configurations

We started our study with the pure dipole case with dimensionless B_⋆^d = 89.3 (recall that R_⋆ = 1). The left panel of Fig. 2 shows the density and the streamlines of angular momentum flux (solid and dashed blue lines) carried by the magnetic field lines at a time of 2T₀, with T₀ the orbital period at R_⋆. The matter accretes through two funnel streams, forming arc-like structures near the magnetic poles. The formation of the funnels occurs at the truncation radius R_trunc ≃ 3.48R_⋆ predicted by Eq. (18), and no signs of accretion are observed in the disk’s midplane. Additionally, the disk undergoes vertical compression from R = R_trunc to R ≃ 10R_⋆, whereas it undergoes a slight thickening beyond 10R_⋆.

Fig. 2. Density maps (in units of ρ0) at t = 2T0 for a dipolar (D model) magnetic configuration (left panel) and for a quadrupolar (Q model) magnetic configuration (right panel). The blue lines on the left side of each panel show the angular momentum flux (fB) carried by the field (see Eq. (16)), with solid blue indicating positive flux and dashed blue indicating negative flux. In the right side of each panel, the yellow vectors show the velocity and the white lines represent the poloidal magnetic field lines. Lastly, inserts show a zoomed-in view of a region of the density map where R < Rco where the gas clearly follows different paths toward the star for the D and Q models. However, it is worth noting that, in both cases, the gas flow toward the star is symmetric about the midplane.

The right panel of Fig. 2 shows the case of a quadrupolar magnetic configuration, which is slightly more complex than the previous one. In this configuration, the streamlines of angular momentum flux carried by the magnetic field lines form mirror lobes relative to the disk’s midplane (see the second panel of Fig. 1). As a result, the angular momentum flux carried by the magnetic field lines converges toward the star, very close to the disk’s midplane. This causes a strong compression of the disk gas for R ≤ R_co and expansion for R_co < R ≤ 10R_⋆. Therefore, at 2T₀, we find a cone-shaped structure near the star that produces gas accretion only in the midplane of the disk. Furthermore, in the coronal region, we find that this magnetic field configuration produces the formation of prolate-spheroid ejections slightly inclined with respect to the midplane.

In both models (D and Q), at t = 2T₀, the magnetic field lines (solid white lines) are modified, and as a consequence, the magnetic field alters the equilibrium state. However, the forces are not strong enough to disrupt the disk.

On the other hand, from Fig. 2 it can be seen that the shape and behavior of the magnetic field lines differ significantly between the dipolar and quadrupolar models. For instance, in the dipolar model, a closed magnetosphere forms near the star with a set of inflated field lines at larger distances. The field lines are observed to inflate and enter through the poles of the star, producing a slight asymmetry in the region of the flattened disk. In contrast, in the quadrupolar case, the poloidal lines enter through the midplane preserving a symmetry in the disk. Finally, we note that in our simulations, inflation of the poloidal lines in the quadrupolar configuration seems somewhat suppressed which could be explained by a relatively dense corona.

3.2. Dipolar-quadrupolar combinations

To model the corona of FS CMa stars, we considered various configurations formed by combining dipolar and quadrupolar fields, as indicated in Table 2 with the letters DQ. The subscript under the logical name “DQ” represents the ratio B_⋆^d/B_⋆^q.

Figure 3 shows the models DQ₁₇, DQ₁₂, DQ₁ and DQ₃, at t = 3T₀. An inspection of this figure reveals that the last three models exhibit a flattened region of the disk that extends from the truncation radius to approximately R ≃ 10R_⋆, as the dipole strength increases. Furthermore, accretion occurs mainly through the formation of a funnel flow below the midplane of the disk (θ > π/2). When the quadrupole polar strength dominates, the flattened region of the disk is considerably reduced and twisted, resembling an extension of the funnel flow. Remarkably, the second funnel flow does not form above the midplane, because the magnetic field lines are more inflated and enter the star closer to the north pole (see the solid white lines in Figure 3 for z > 0).

Fig. 3. Gas density maps for DQ models (see Table 2) at t = 3T0. We overlapped the yellow vectors to show the direction of the gas velocity and the poloidal magnetic field lines in each case (solid white lines). Note that the gas follows the magnetic field lines forming only a twisted funnel flow that enters the star below the midplane of the disk in all cases. Here the inserts show a zoomed-in view of a region of the density maps contained within the corotation radius.

In the particular model DQ₁₇, we find that, as the quadrupole magnetic polar strength dominates over the dipole, the gas tends to follow the magnetic field lines closer to the disk’s midplane. However, the effect of the dipole is still visible, as a twisted funnel flow forms for R < 10R_⋆ above the disk’s midplane, leading to an inflated region at the edge of the disk. Note that although this funnel forms above the disk, it enters the star below the midplane (see the insert in the first panel of Fig. 3). In fact, in this model we can see the formation of two gas lobes, one above and one more elongated below the midplane. Additionally, unlike in the pure dipole case, here the gas flow in the funnel can vary because the path traced by the field lines changes due to their expansion, resulting in variable accretion of matter (see Section 4.2.1). In other words, in this dipole-quadrupole combination, intermittent gas accretion can occur; this latter accretion could be compared with observations. Finally, the substructures that form in the corona region, reminiscent of wavy magnetospheric ejections, also become more asymmetric as the quadrupole polar strength dominates.

3.3. Complex magnetic configuration

We conducted simulations using a dipolar plus octupolar magnetic field, labeled as DO₁₇ and DO₁₃ in Table 2. Let us first analyze the case DO₁₃.

Figure 4 shows the temporal evolution of the gas density. At t = 3T₀, a pair of well-defined funnels with a crab-claw topology forms, channeling material onto the star near the disk midplane. Additionally, two much fainter funnels are observed near the stellar poles. By t = 6T₀, multiple new funnels have developed, stacked within the original crab-claw structures, resulting in a total of six distinct funnels directing material onto the star at different latitudes. The latter result is a consequence of the topology of the magnetic field lines in the octupole (see the third panel of Fig. 1). We stress that the funnels at the top are thinner and more tenuous. The density contrast between the top and bottom of the funnels is also reflected in the coronal region, as dense asymmetrical arches form below the midplane. Finally, we note that the lower outer funnel resembles a dipole shape.

Fig. 4. Temporal evolution of the gas density (in logarithmic scale in units of ρ0) for the model DO13. Left panel : Time t = 3T0. Right panel: Time t = 6T0. The streamlines of the angular momentum fluxes (fB) carried by the magnetic field are represented by the blue lines in the left side of each panel. On the right sides the velocity vectors are shown with yellow arrows. The inserts show a zoomed-in view of a region of the density maps where R < Rco.

On the other hand, for model DO₁₇, we observe again the formation of multiple funnels following a pattern very similar to the magnetic field lines shown in the lower panel of Fig. 1. In this case, the magnetic field strength of the octupolar component is seven times greater than that of the dipolar component. As a result, the accretion funnels above and below the midplane of the disk are both well defined. Nevertheless, the inclusion of the dipolar component introduces a twist in the region where each funnel originates, ultimately leading to an asymmetry between the upper and lower regions of the disk and corona (see the left panel of Fig. 5).

Fig. 5. Gas density comparison between dipolar plus octupolar models in logarithmic scale in units of ρ0 at t = 5T0. Left panel: Model DO17. Right panel: Model DO13. The blue lines and inserts are the same as in Fig. 4. The gas flow is compared at t = 5T0, when a quasi-steady state has been reached in the mixed dipolar–octupolar configurations.

In Fig. 5 we compared models DO₁₃ and DO₁₇ at t = 5T₀. Far from the star (R > 10R_⋆), the disk structure is similar in the two cases. However, for R ≤ 10R_⋆ and in the regions where the funnels form, there is clearly a substantial difference. As the strength of the dipolar magnetic field increases (model DO₁₃), the disk is compressed. Interestingly, the gas accumulates to a greater extent below the midplane, which may lead to a different accretion rate compared to the quasi-octupolar case (DO₁₇). In both models, the perturbations in the corona are asymmetric with respect to the midplane (an effect that can also be inferred from the streamlines of angular momentum; see the blue lines in Fig. 5), with the model DO₁₇ exhibiting a stronger perturbation.

4. Discussion

4.1. Observational motivation for the different magnetic configurations

It has been suggested that the magnetic configuration on the surface of the stars may be a composition of dipolar, quadrupolar and even octupolar configurations (e.g., Zanni & Ferreira 2009; Čemeljić 2019). On the other hand, observations suggest that the magnetic field on the surface of stars may be more complex than the dipolar configuration usually considered in numerical studies. Preliminary spectropolarimetric results for the observations of one FS CMa star, IRAS 17449+2320, show a strong magnetic field with a predominant dipolar configuration plus quadrupolar contribution. The estimated dipolar field strength is around 9000 G, while the quadrupolar component is around 3000 G (Bermejo-Lozano et al., in prep). For this reason, to model the corona of FS CMa-type stars, it is reasonable to assume different configurations of the magnetic field on the stellar surface.

4.2. Twisted funnel flow accretion

As shown in Section 3, the accretion onto the star depends on the relative contributions of the multipolar components. Depending on the dominant magnetic polar strength in the configuration, accretion can occur through two funnel flows that are completely symmetric with respect to the disk midplane (purely dipolar configuration), or it can be drastically modified (see Section 4.2.3 for further discussion).

In a DQ configuration, regardless of the B_⋆^d/B_⋆^q ratio, a single twisted funnel flow is formed, connecting the star below the disk midplane. That is, gas accretion onto the star occurs below the stellar equator. It is important to emphasize that the twisting of this single funnel flow becomes more pronounced as the quadrupole magnetic moment increases. Under this scenario, gas flows are generated above the disk midplane, which fall toward the stellar surface (see Fig. 3). It is inferred that a twisted funnel flow can significantly change the mass flux of gas entering the star. Note that in both the twisted funnel flow and the other cases studied here, the quasi-steady state is reached within a few star periods (t ≤ 5T_⋆). With this in mind, in the next subsection we analyze the accretion of gas onto the star.

4.2.1. Mass accreted in a disk

The mass flow rate crossing the whole polar range is:

$$\begin{array}{c}\hfill \dot{M}=-4\pi R^2{\displaystyle {\int }_0^{\pi /2}}\rho v_{R}sin\theta \mathrm{d}\theta .\end{array}$$(20)

Figure 6 shows the mass flow through the disk onto the surface of the star R = R_⋆ and at the truncation radius R_trunc = 3.48R_⋆ given by Eq. (18), for the models DQ₁₇, DQ₁₂, DQ₁ and DQ₃. In all four models, the magnitude of mass flux at R = R_⋆ increases with the strength of dipolar component. In other words, when the intensity of the dipolar magnetic field increases, Ṁ also increases until it reaches a quasi-stationary value (which occurs in about 5 orbital periods of the star). From this figure, we can infer the dependence of Ṁ on the path the gas follows toward the star, which is constructed based on the initial magnetic configuration. For instance, for the models DQ₁ and DQ₃, Ṁ takes a very similar value (∼0.1) at t > 4.5T₀, which is, in fact, higher than it is for the other two models (where Ṁ ≃ 10⁻³), with the quadrupole component making a larger contribution. The above argument can be confirmed from the funnel structures shown in Fig. 3. The funnels in models DQ₁ and DQ₃ have a thinner structure, and the footpoint of these funnels is closer to the midplane of the disk. This feature is a consequence of the strong compression of the disk between R ≃ R_trunc and R ≃ 10R_⋆. Therefore, the distance traveled by each fluid element is shorter in these cases than in models where the funnel arch is increased, such as DQ₁₇ and DQ₁₂.

Fig. 6. Time dependence of the mass flow rate for DQ models at R = R⋆ (upper panel) and at the truncation radius R = 3.48R⋆ (lower panel).

Consider now the flux through the sphere with radius R = 3.48R_⋆ (that is, the truncation radius given by Eq. (18) for a pure dipolar configuration). A faster flow toward the star occurs in the DQ₁₇ configuration within the first three orbital periods. This is expected because, although a twisted funnel forms, the dominance of the quadrupolar magnetic component leads to a greater gas flow above and below the disk midplane (see the first panel of Fig. 3). Nevertheless, for t ≥ 3T₀, the mass flow decreases in a stepwise manner, corresponding to episodic expansions in the length of the funnel flow (due to the expansion of the magnetic field lines). Something similar happens in the model DQ₁₂ when t ≥ 2T₀.

On the other hand, when the dipolar magnetic moment is equal to or greater than the quadrupolar magnetic moment, as in the models DQ₁ and DQ₃, the mass flux toward the star remains quasi-stationary for t > 1T₀, which means that the funnel is more stable and does not undergo strong radial expansion. The behavior in the model DQ₁ (where the intensity of the dipolar and quadrupolar magnetic fields are equal) suggests that the angular momentum flux is mainly carried by the magnetic field. To verify this hypothesis, we considered the angular momentum fluxes carried by the field f_B, and by the matter f_m (see Eqs. (16) and (17)). Figure 7 shows the maps of f_B and f_m for model DQ₁ at t = 4T₀. It can be seen that the angular momentum flux driven by the field is more negative in the region where the funnel formation takes place (see the left panel) compared to the angular momentum flux f_m driven by the matter (see the right panel). This produces a radial shift of the funnel flow that inevitably results in a sharp drop in the mass flow Ṁ across the truncation radius (see the lower panel of Figure 6).

Fig. 7. Maps of the angular momentum fluxes carried by the field (fB) and the matter (fm; see Eqs. (16) and (17)). Note that the arc formed by the matter flux (fm; right) is contained within the region of negative angular momentum flux (fB) through the magnetic field (left).

4.2.2. Tracking the quadrupolar truncation radius

As specified in Section 2.3.2, the theoretical value of the truncation radius for the purely dipolar and octupolar cases is determined by Eqs. (18) and (19), respectively. However, due to the nature of the initial configuration of the magnetic field components, which implicitly depend on a poloidal configuration, it is not possible to analytically predict the value of the truncation radius for the quadrupolar case at the disk midplane.

To obtain an expression for R_trunc for the quasi-quadrupolar case, we calculated the truncation radius from the simulation results for each of the DQ models and inserted it into function of the intensity of the dipolar component of the magnetic field. From these data we obtained the next fitting formula:

$$\begin{array}{c}\hfill \frac{R_{\mathrm{trunc}}}{R_{\star }}=K_D{\left(B_{\star }^d\right)}^2-K_Q{B}_{\star }^d+K_0\end{array}$$(21)

where K_D = 2.37 × 10⁻⁴, K_Q = 3.78 × 10⁻² and K₀ = 4.97.

Figure 8 shows the comparison of the fitting formula given in Eq. (21) and the simulation data. For reference, we have included the value of the purely dipolar case predicted by Eq. (18). We see that our formula adequately predicts the truncation radius when the quadrupolar component is considerably dominating.

Fig. 8. Truncation radius at t = 3T0 for DQ models. The intersection point of the dashed purple and solid orange lines represents the theoretical truncation radius value (R = 3.48R⋆) calculated using Eq. (18).

4.2.3. Beyond the dipolar funnel flow

In the standard funnel flow formation, the gas falling onto the star follows the path defined by the magnetic field lines. In general, we find this behavior across all the magnetic configurations studied. However, the fact that the gas follows the magnetic field does not necessarily imply that a funnel flow is formed. For instance, a purely quadrupolar configuration produces an accretion flow confined to the disk midplane, thereby preventing the lifting of material above the disk midplane and formation of a funnel flow (see the right panel of Fig. 2). Furthermore, we find multiple funnel formation for octupolar cases.

Note that even for a longer time evolution (which is beyond the scope of this work) in both the octupolar configurations and the other magnetic configurations analyzed in this study, the inflation of the magnetic field lines is supported by the differential rotation between the star and the gas, as well as the value of the magnetic diffusivity in the disk (extended magnetosphere; see Moranchel-Basurto et al. 2024, and references therein for details). Furthermore, in a dipolar magnetic configuration including background magnetic field, if the constraint of a thin disk is relaxed, accretion bursts, i.e. gas streaming directly from the truncation radius to the star induced by the interchange instability, can be found as well as one-sided outflows for stars with kilogauss magnetic fields. In such cases magnetospheric accretion develops (see Gaches et al. 2024, for details). Therefore, the gas accretion rate and density structures (including the twisted funnels formed in DQ models) can remain stable longer than reported here.

5. Conclusions

We performed 2.5-D MHD simulations of star-disk magnetospheric interaction with the aim of determining the gas density structure and the gas accretion process in massive stars with strong magnetic fields, a characteristic of several Herbig Ae/Be and FS CMa stars.

We have found that the gas accretion process is highly dependent on the symmetric (asymmetric) configurations of the magnetic field. We can classify accretion flows as follows:

• (i)

  In a dipolar (symmetric) magnetic configuration, two funnel streams that are symmetric with respect to the disk’s midplane form.

• (ii)

  In a quadrupolar (symmetric) magnetic configuration, accretion is driven along the disk’s midplane, which is flattened by the quadrupolar magnetic field lines. This process is accompanied by a widening of the disk due to the conical shape of the stream, with the base of the cone located near the truncation radius.

• (iii)

  In the case of a dipole plus quadrupole (asymmetric with respect to the disk midplane) magnetic configuration, we find the formation of twisted funnel flows when the quadrupolar polar strength component dominates. These funnel streams connect with the surface of the star below the disk midplane.

• (iv)

  Finally, for the dipole plus octupole configuration, we find the formation of multiple funnels behind the truncation radius predicted by Eq. (19). Remarkably, the inclusion of the dipolar component produces an asymmetry in the funnel structure with respect to the disk midplane.

In addition, in all our models, we find density substructures in the corona that follow different patterns depending on the magnetic field configuration, with these structures being more pronounced in multipolar configurations. Since a multipolar magnetic field configuration can be the result of a post-merger event, this suggests that, for a strong magnetic field such as in the case of several FS CMa stars, the disk and corona can exhibit a complex density distribution as consequence of a post-merger event. In future work, we plan to investigate the cases where magnetic field configurations at the surface of the star are misaligned with respect to the rotational axis (Ω). To do so, full three-dimensional models will be considered.

Acknowledgments

We thank the referee for useful suggestions. The work of R.O.C. was supported by the Czech Science Foundation (grant 21-23067M). Computational resources were available thanks to the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). MČ acknowledges the Czech Science Foundation (GAČR) grant No. 21-06825X and the Polish NCN grant 2019/33/B/ST9/01564.

References

Appendix A: Vector potentials and flux functions

We assume that the star is a uniformly magnetized rotating sphere, on which flow the azimuthal currents I(θ). The distribution of currents I(θ) will depend on the magnetic configuration so that if it is a purely dipolar magnetic field I(θ) = I_d sin θ, in the case of a purely quadrupole magnetic field I(θ) = I_q sin θ cos θ, while when it is a purely octupolar configuration we have I(θ) = I_oc sin θ (5cos²θ − 1) (Shakura et al. 1991). Since we use a simplified axisymetric model, the vector potential has only azimuthal component. Therefore, the potential vector for dipolar, quadrupolar and octupolar configurations reads

$$\begin{array}{c}\hfill A_{\varphi }^d\equiv \frac{{\mu }_d}{R^2}sin\theta ,\end{array}$$(A.1)

$$\begin{array}{c}\hfill A_{\varphi }^q\equiv \frac{3{\mu }_q}{4R^3}sin\theta cos\theta ,\end{array}$$(A.2)

$$\begin{array}{c}\hfill A_{\varphi }^{\mathit{oc}}\equiv \frac{{\mu }_{\mathit{oc}}}{2R^4}sin\theta (5{cos}^2\theta -1).\end{array}$$(A.3)

By definition, B = ∇ × A and, after some algebra, we can derive the components of the magnetic field.

The contours of the magnetic flux function

$$\begin{array}{c}\hfill \psi =A_{\varphi }Rsin\theta \end{array}$$(A.4)

trace the poloidal field lines. The magnetic flux functions for the dipolar, quadrupolar and octupolar fields are given by

$$\begin{array}{c}\hfill {\psi }^d=\frac{{\mu }_d{sin}^2\theta }{R},\end{array}$$(A.5)

$$\begin{array}{c}\hfill {\psi }^q=\frac{3{\mu }_q}{4R^2}{sin}^2\theta cos\theta ,\end{array}$$(A.6)

and

$$\begin{array}{c}\hfill {\psi }^{\mathit{oc}}=\frac{{\mu }_{\mathit{oc}}{sin}^2\theta (5{cos}^2\theta -1)}{2R^3}.\end{array}$$(A.7)

All Tables

All Figures

	Fig. 1. Initial density (in units of ρ0) and magnetic field lines (white curves) for model D (first panel), model Q (second), model O (third), model DQ (fourth) and model DO (fifth). In all cases, the sum of the different multipolar contributions to the stellar magnetic field is equal to 89.3. Here, r and z represent the cylindrical coordinates, measured in units of the stellar radius.
In the text

Fig. 2. Density maps (in units of ρ0) at t = 2T0 for a dipolar (D model) magnetic configuration (left panel) and for a quadrupolar (Q model) magnetic configuration (right panel). The blue lines on the left side of each panel show the angular momentum flux (fB) carried by the field (see Eq. (16)), with solid blue indicating positive flux and dashed blue indicating negative flux. In the right side of each panel, the yellow vectors show the velocity and the white lines represent the poloidal magnetic field lines. Lastly, inserts show a zoomed-in view of a region of the density map where R < Rco where the gas clearly follows different paths toward the star for the D and Q models. However, it is worth noting that, in both cases, the gas flow toward the star is symmetric about the midplane.

In the text

Fig. 3. Gas density maps for DQ models (see Table 2) at t = 3T0. We overlapped the yellow vectors to show the direction of the gas velocity and the poloidal magnetic field lines in each case (solid white lines). Note that the gas follows the magnetic field lines forming only a twisted funnel flow that enters the star below the midplane of the disk in all cases. Here the inserts show a zoomed-in view of a region of the density maps contained within the corotation radius.

In the text

Fig. 4. Temporal evolution of the gas density (in logarithmic scale in units of ρ0) for the model DO13. Left panel : Time t = 3T0. Right panel: Time t = 6T0. The streamlines of the angular momentum fluxes (fB) carried by the magnetic field are represented by the blue lines in the left side of each panel. On the right sides the velocity vectors are shown with yellow arrows. The inserts show a zoomed-in view of a region of the density maps where R < Rco.

In the text

	Fig. 5. Gas density comparison between dipolar plus octupolar models in logarithmic scale in units of ρ0 at t = 5T0. Left panel: Model DO17. Right panel: Model DO13. The blue lines and inserts are the same as in Fig. 4. The gas flow is compared at t = 5T0, when a quasi-steady state has been reached in the mixed dipolar–octupolar configurations.
In the text

	Fig. 6. Time dependence of the mass flow rate for DQ models at R = R⋆ (upper panel) and at the truncation radius R = 3.48R⋆ (lower panel).
In the text

	Fig. 7. Maps of the angular momentum fluxes carried by the field (fB) and the matter (fm; see Eqs. (16) and (17)). Note that the arc formed by the matter flux (fm; right) is contained within the region of negative angular momentum flux (fB) through the magnetic field (left).
In the text

	Fig. 8. Truncation radius at t = 3T0 for DQ models. The intersection point of the dashed purple and solid orange lines represents the theoretical truncation radius value (R = 3.48R⋆) calculated using Eq. (18).
In the text