2.5-MHD models of circumstellar discs around FS CMa post-mergers : II. Stationary accretion stage

We study the star-disc interaction in the presence of the strong magnetic field ( 𝐵 ★ = 6 . 2 𝑘𝐺 ) of a slowly rotating star. This situation describes a post-merger of the spectral type B and has not been previously investigated. We perform a set of resistive and viscosity 2 . 5 𝐷 -magnetohydrodynamical simulations using the PLUTO code. Based on our previous work, we consider the initial gas disc density 𝜌 𝑑 0 = 10 − 13 gcm − 3 since it describes the conditions around IRAS 17449+2320 well. We find that the fall of gas towards the star occurs in the mid-plane, and remarkably, intermittent backflow takes place in the mid-plane in all of our models for 𝑅 ≥ 10 𝑅 ★ . However, we do not rule out that the funnel effect may occur and cause the accretion closer to the poles. Also, when larger values of viscosity ( 𝛼 𝜈 = 1) and stellar rotation rate ( 𝛿 ★ = 0 . 2) are considered, we find that the disc exhibits a thickening which is characteristic of FS CMa-type stellar objects. Additionally, we find that the poloidal magnetic field lines twist over short periods of time, leading to magnetic reconnection causing coronal heating that could explain the presence of the Raman lines found observationally in several FS CMa stars. Lastly, we find the formation of several knots in the magnetic field lines near and in the mid-plane of the disc which produce perturbations in the density and velocity components, as well as the formation of shallow gaps whose position depends on the inflation of the magnetic field lines.


INTRODUCTION
In the recent two decades the theoretical and numerical effort to understand the dynamics of gas during the interaction between the magnetosphere and the inner region of an accretion disc of a rotating-magnetized star has produced important advances for some types of stars such as: classical T Tauri stars, white dwarfs, neutron stars, X-ray pulsars, among others (Koldoba et al. 2002;Romanova et al. 2002;Long et al. 2005;Bessolaz et al. 2008;Kluźniak & Rappaport 2007;Zanni & Ferreira 2009, 2013;Čemeljić 2019;Ireland et al. 2021;Čemeljić et al. 2023).Among the most important results are the bounding of the disc truncation radius value due to the star's magnetosphere and the formation of accretion funnel flows.The latter is implicitly associated with the spin-up torque felt by the star and represents a stationary gas accretion state.However, both outcomes remain in the regime of sun-like mass stars including a weak magnetic field ( ★ ≤ 1).
Less studied is the star magnetosphere-disc interaction regime for more massive stars with a strong magnetic field (Miller & Stone 1997;Moranchel-Basurto et al. 2023;Čemeljić et al. 2023).Apart from the expensive computational cost, the study of this type of interaction represents a complex problem itself, as has been well pointed out in Moranchel-Basurto et al. (2023, hereinafter called Paper I), since very strong magnetic field coming from the central star can produce a highly dynamic accretion of gas (non-stationary stage).
★ E-mail: amoranchel087@gmail.comIt is worth mentioning that the stationary or non-stationary regime of the accretion depends mostly on the initial conditions and the strength of the magnetic field does not play the crucial role.Therefore, being able to restrict the free parameters to the smallest possible number can make the study of the interaction of the magnetosphere with the inner disc more tractable, at least in numerical models.For instance, Čemeljić (2019) studied star-disc magnetospheric interaction (including only a high viscosity in the disc,   = 1) by an 'Atlas' of numerical models trying to obtain quasi-stationary solutions which can be applied to accretion discs around of the YSO-type objects with a low magnetic field ( ★ < 1).He found that increasing stellar magnetic field in the simulation, the gas density of the disc increases.Furthermore, he found that in models with a faster rotation of the star, additional solutions are produced: (i) a fast axial outflow, and, (ii) with both the conical and axial outflow, similar to what was reported in Romanova et al. (2009).With all the above in mind, the parameters obtained from the observations play an important role.Within the regime of massive stars with a strong magnetic field are the objects called FS CMa post-mergers.Particular attention must be paid to the member of this set labeled as IRAS 17449+2320 object.Because in a recent study the magnitude of the magnetic field on the surface of the star was determined ( ★ = 6.2; see Korčáková et al. 2022, and Paper I for details).So, similar to Paper I, we are interested in studying the star-disc magnetosphere interaction in said regime taking into account the observational parameters for the FS CMa stars, such as: the mass ( ★ = 6 ⊙ ) and stellar radius ( ★ = 3 ⊙ ) (Kříček et al. 2017), and the value of the magnetic field mentioned above.Nevertheless, in this case with the aim of seeing if the stationary accretion stage takes place, unlike Paper I, here we include as free parameters the stellar rotation rate  ★ , the magnetic diffusivity  m and the viscosity   within the disc in our numerical models.
The paper is laid out as follows.In Section 2 we present the physical model, while the numerical implementations used in our 2.5D magnetohydrodynamical simulations are discussed in section 3.In Section 4 we present the results of our numerical models.A brief discussion is presented in Section 5 and finally concluding remarks can be found in Section 6.

PHYSICAL MODEL
Our physical model is composed of a slightly sub-Keplerian disc around a magnetized rotating star with a dipolar magnetic field configuration.The models presented in this study are numerical solutions of the magnetohydrodynamical (MHD) equations in which the effects of viscosity and resistivity within the disc are taken into account.Apart from the free parameters considered in this study, our physical model follows those presented in Paper I.

MHD equations
The magnetohydrodynamic equations governing the gas dynamics are the conservation of mass the conservation of momentum equation

𝜕 𝜌v 𝜕𝑡
the energy equation the last term on the left side of the energy equation (3) includes the heating terms, to balance this, note that the cooling term Λ  has been included on the right side of the equation, that means, that the system should evolve adiabatically.Said cooling term also helps prevent the thermal thickening of the gas disc (Zanni & Ferreira 2009).
And finally the induction equation by which the description of the evolution of the magnetic field is described is given by: where  is the thermal gas pressure,  the gas density, v its velocity, I the unit tensor, B the magnetic flux density vector, J is the electric current density given by the Ampére's law J = ∇ × B/4,   is the resistivity and is related with the magnetic diffusivity by   =   /4, Φ =   * / the gravitational potential and  represents the total energy density given by where  = 5/3 is the polytropic index of the plasma.On the other hand, the viscous stress tensor  is defined as: for the BSROT 2b model with a higher rate of stellar rotation  ★ = 0.2.In this case we find the formation of multiple disc layers between  = 5 0 and  = 10 0 as well as the inflation of the disc that resembles the formation of a huge funnel.The white lines show the poloidal magnetic field lines.Note that in both models except for the increase in the stellar rotation rate in a factor 2, all other parameters remained fixed.
where   is the dynamic and   =   / the kinematic viscosity, respectively.

Accretion disc
The density   and pressure   of the gas in the accretion disc is set following the three-dimensional models of Keplerian accretion discs considering spherical coordinaties (, , ) used in Zanni & Ferreira (2009) (which in turn is based on the disc model of Kluzniak & Kita 2000): where ℎ =   /  is the aspect ratio, with   and   = √︁   * / the sound speed and the Keplerian velocity, respectively.In Eq. ( 7)  0 and  0 denote the initial values of the density and the Keplerian velocity calculated at  =  0 in the mid-plane of the disc.Note that we define  =  sin  as the cylindrical radius.

Disc atmosphere
In a similar form to Zanni & Ferreira (2009), we included a non-rotating polytropic hydrostatic atmosphere with a density with  = 5/3.The density contrast between the disc and the atmosphere is  0 atm / 0 = 0.01, which is kept fixed in all our models.

Magnetic field configuration
We assume that the star is located at the origin of the coordinate system.The stellar magnetosphere is modeled initially as a purely dipolar field aligned with the rotation axis of the star-disc system.The magnetic field is defined by the flux function where  * is the magnetic field at  * and  = 0.The radial and polar field components are therefore given respectively by and The relation between flux function and potencial vector is given by: and then the components of potential vector in spherical coordinates are given by: = 0 and   = 0.It should be noted that for a suitable numerical treatment of the magnetic field, we employ the "field-splitting" technique (Tanaka 1994;Powell et al. 1999) which is widely used in this type of studies (see Čemeljić 2019, and references therein).

Magnetic diffusivity models
The role of magnetic difussivity   can be important, this can arise from turbulence generated in the disc and this parameter may be responsible for modifying the strength of the magnetic field lines.For this reason we decided to explore different descriptions of resistivity   = 4  , considering that the magnetic diffusivity   is given by: • Classical description (Zanni & Ferreira 2009;Čemeljić 2019): • Modified first version (see Masoumeh & Khesali 2022, and references therein) • Radial and vertical dependence as given in Bessolaz et al. (2008): the latter decreasing on a disc scale height  =   /Ω  , here Ω  is the Keplerian rotation speed.In Eq. ( 20),   is a constant that we fixed to the value 0.1.Throughout this study we will focus in more detail on the resistivity given in Eq. ( 21), since this expression takes into account the vertical stratification of the accretion disc.We have taken the -axis parallel to the rotation axis of the disc and assume that the magnetic dipole moment is aligned with the rotational axis.

Code
In order to investigate the disc-star interaction, we use the MHD module provided by PLUTO code version 4.4 (Mignone 2009, Mignone et al. 2007) to perform adiabatic global MHD simulations of accretion discs to numerically solve the MHD equations given in section 2.1.Simulations runs on spherical geometry (, , ), were performed using the second-order linear reconstruction Runge Kutta integration in time to advance conserved variables.To enhance the stability of the code, in the subroutine pml states we consider the Van Leer limiter in the density and in the magnetic field since is more diffusive, while a monotonized central differences limiter (minmod) was used in pressure and velocity fields.Also, we mention that for the inclusion of the cooling term that controls the possible heating of the disc, we have used a power-law cooling option provided in the code.Additionally, inspired by the work of Čemeljić (2019) we use a hll solver with a modification in the flag shock subroutine: flags were set to switch to a more diffusive hll solver when the internal energy defined by  =  gas /( − 1) was lower than 1% of the total energy given by the equation ( 5).Furthermore the condition ∇ •  = 0 was maintained using the constrained transport method.Lastly, we mention that in this work we have followed closely the initial conditions in the disc and corona and boundary conditions at the edges of the computational domain as in the works of Čemeljić (2019) and Zanni & Ferreira (2009), in the next subsections we described these conditions in more detail.

Simulation parameters and units
The parameters of the simulations are the same as were used in the numerical models of the paper I, and were selected to match a set of observed values for the FS CMa stars (Korčáková et al. 2022).More specifically, we fix the mass, radius, initial disc density and the magnetic field of the star to be  * = 6 ⊙ ,  * = 3 ⊙ ,  = 1 × 10 −13 g cm −3 and  * = 6.2 × 10 3 .Furthermore we take into account that the aspect ratio is ℎ = 0.1.The different models considered are presented in Table 1.
Once the initial conditions are normalized, the problem depends on six dimensionless parameters: the aspect radio ℎ, the equatorial stellar field intensity  * / 0 , the stellar rotation rate  * =  * Ω * /  * , the coronal density contrast  atm / 0 , and the viscous and resistive coefficients   and   .Except for the transport coefficients and the initial differential rotation of the star, the other parameters are the same as were used in the simulations of paper I:

Co-rotation and truncation radius
To study the interaction of a star-magnetized disc it is important to identify where the truncation radius must occur by the stellar magnetosphere.So far there are numerous studies that debate the most precise way to determine this radius (see for instance : Romanova et al. 2002;Zanni & Ferreira 2009;Bessolaz et al. 2008;Long et al. 2005).However, despite the differences between all the models, it has been determined that for accretion to occur the truncation radius   must be located within the co-rotation radius  co = (  ★ /Ω ★ ) 1/3 .We can rewrite the co-rotation radius in terms of the star radius  ★ if we isolate Ω ★ variable from the definition of stellar rotation rate  ★ given in the previous section and substitute its value in the definition mentioned above, we will obtain after doing a bit of algebra that the radius of co-rotation can be written as follows  co = 4.64 ★ .
On the other hand, we obtained the value for the truncation radius following the next expression (Bessolaz et al. 2008): considering that the value of   can be close to 1 and substituting the values described in section (3.2) we obtain that: It should be emphasized that in the calculation of   given in Eq. ( 23), we have considered an accretion rate  = 10 −8  ⊙ yr −1 , which is a conservative value because there is currently debate about what is the correct value of the accretion rate in the FS CMa-type objects (see Carciofi et al. 2010, and references therein).

Mesh domain
The axisymmetric disc is modeled in the radial range 1 0 ≤  ≤ 30 0 with a polar angular extension in one quadrant (0 ≤  ≤  2 ).All of our numerical simulations presented here were performed by the same numerical resolution, with a logarithmically increasing grid in radial direction   = 109 and,   = 50 polar grid cells uniformly distributed.The resulting computational mesh is sufficient to resolve properly the star-disc magnetospheric interaction (SDMI) in both dimensions (see Čemeljić et al. 2017), and in turn allows us to study a larger number of parameters (see Table 1) with a reasonable computational cost.

Boundary Conditions and internal boundary
Since we are interested in analyzing the gas dynamics when SDMI occurs, considering a slow-rotating massive star with strong magnetic field (see subsection 3.2), we take particular care about the boundary conditions implemented on the stellar surface as in previous studies where the SDMI has been analyzed for young stars with a weaker magnetic field (see for instance Romanova et al. 2002Romanova et al. , 2009;;Zanni & Ferreira 2009;Čemeljić 2019).Then, the boundary conditions that were used in carrying out this work closely follow those described in the works of Zanni & Ferreira (2009) and Čemeljić (2019).
To obtain adequate gas dynamics on the surface of the star, we implement the following conditions.First, in the reference frame co-rotating with the star, the electric field is zero (see equation (11) in Zanni & Ferreira 2009, for details).Second, to prevent low-density near the star, we modified the density value when it falls below some limit value 1 within a grid cell of the computational domain on the surface of the star.For the corona region, to conserve the same sound speed and to keep the conservation of the momentum we modify the gas pressure and the velocities, respectively.Lastly, to follow the evolution of the disc, we activate the scalar tracer, taking a value of unity in the disc domain and zero outside of it.Note that because the resistivity described in the equations ( 19)-( 21) is only applied within the domain of the disc mesh, the reconnection of the magnetic field arising in the magnetosphere is handled by numerical dissipation.
For the inner grid ghost cells, following Čemeljić (2019) we prescribe the values of the density, pressure, and toroidal components of the velocity and magnetic field from the active zones.This implies that the density and the magnetic field are linearly extrapolated considering a Van Leer limiter, and minmod limiter in the pressure and velocity (see subsection 3.1).It should be emphasized that to achieve numerical stability in the corona, when   > 0 in the inner radial boundary, we define the gas pressure as  = (2/5−   2 )    −5/2 , with   a free parameter (see Čemeljić 2019).
On the other hand, in the radial outer boundary condition, we consider the outflow condition for the velocity components, and for the density and magnetic field a linear extrapolation was used.Lastly, in the -direction we use the PLUTO code options 'axisymmetric' and 'eqtsymmetric' for the boundary conditions at  =  min and  =  max , respectively, since we simulate only one hemisphere of the disc.

RESULTS
In this section we present the results of our 2.5-MHD simulations of a resistive accretion disc for different values of   and  m , and also considering different stellar rotation rates.Our study is focused on the strong magnetic field regime measured in massive stellar objects such as FS CMa-type stars (see Korčáková et al. 2022).
All the parameters and the different models are given in the table 1.The viscosity parameter describes the strenght of the viscous torque, that allows the disc to accrete, while the resistivity parameter   allows coupling of the stellar magnetic field with the disc material.
For both viscosity and resistivity we consider the pair of values of 0.1 and 1.So the Prandtl number can take the set of values of  m = {0.06,0.6, 6.6}.The motivation for choosing these values in the viscosity lies in the fact that previous studies of the star-disc interaction show that if it is fulfilled  <  crit ≃ 0.685, a backflow is originated in the mid-plane of the disc, which results in the accretion of the disc towards its surface (see Zanni & Ferreira 2009, and references therein).
On the other hand, the choice of values 0.1 and 1 for the  m resistivity coefficient, are motivated to know which region of the disc is steadily connected to the star in our models with a strong magnetic field.Since the opening of the magnetosphere is intrinsically driven by the differential rotation between the disc and the star, as well as with the resulting buildup of toroidal magnetic pressure.The latter can be modified by the value that  m takes on the disc, since if  m = 1 limits the growth of a toroidal magnetic field, producing a larger connected region (an extended magnetosphere).In the case when  m = 0.1 it has the opposite effect on the magnetic pressure (due to the strong coupling), thus producing a smaller connected region (compact magnetosphere).

Comparison of the stellar rotation rate for models with an extended magnetosphere
Figure 1 shows a comparison between the BSROT0 1b and BSROT0 2b models.In these models the only difference in the parameters is a factor of 2 in the rate of stellar rotation  ★ .We find that in the case of a slower stellar rotation (BSROT0 1b model) it produces an accretion towards the star in the mid-plane of the disc as well as a distribution of the gas of the disc outside the mid-plane that resembles the formation of the funnel effect that can give rise to accretion near the poles.Also, we find in this case the formation of disc winds and magnetospheric ejections.In the case of T Tauri stars with a smaller magnetic field on the star's surface ( ★ = 1.1) a similar behavior has been found in the formation of the funnel effect and in the generation of magnetosphere ejections (see Zanni & Ferreira 2013).However, in said study there are no signs of gas accretion in the mid-plane of the disc.In the case of a star rotating a little faster (case described here by the BSROT 2b model), we find an inflation of the disc, which is the result of the inflation and reconnection of the magnetic field lines, in addition to the fact that we find the formation of layers between  = 5 0 and  = 10 0 .
These layers fail to form multiple funnels towards the star due to the increase in the rotation of the star.

Models with compact and extended magnetospheres with the same stellar rotation rate
In the Fig. 2 we show the logarithmic gas density maps for the BSROT0 1a and BSROT0 1b models calculated at  = 60 ★ .In the first case, we find that in a disc with lower magnetic diffusivity ( m = 0.1) there is a formation of two knots at  = 17 and  = 25 0 .
The first knot creates strong disturbances in the density of the disc, producing a partial gap in  = 17 0 .The second knot located outside the mid-plane of the disc (at  = 25 0 ) generates a magnetosphere ejection.For the BSROT0 1b model with a higher resistivity ( m = 1), we find the formation of a single knot outside the mid-plane of the disc at  = 13 0 which produces a magnetosphere ejection.In both models we find that the gas accretion takes place in the mid-plane of the disc.We emphasize that, since the stellar rotation rate in these models is the same, it is clear that the differences in the gas density distribution arise due to the value of the resistivity within the disc.
Figure 3 shows the radial profiles of the sonic Mach number   = |  |/ √︁ /| =0 , the disc angular velocity Ω and the Keplerian velocity Ω K , respectively (nomalized by Ω ★ ), calculated at  = 10 ★ for the BSROT0 1a model.The behaivor of the dynamics in the disc is changing the accretion Mach number which increases and reaches a maximum value of   ≈ 1.3 in a radius of 5 0 (see Fig. 3).Note that   starts to increase outward from the truncation radius   (green vertical line in Fig. 3).On the other hand, we can see in Fig. 3 that the angular velocity of the disc Ω becomes super-Keplerian for  > 3.1 0 (solid blue line).This implicitly means that the rotating star exerts a positive torque on the disc.In other words, the star loses angular momentum which consequently translates into a lower rotation speed of the star.Although in this case, the magnetic field is strongly coupled to the disc (since  m = 0.1), we point out that we find very similar results in cases when the resistivity is higher ( m = 1).Therefore, if the magnetic field is able to couple to the disc, the rotation speed of the star will inevitably decrease.The latter may explain why FS CMa-type stellar objects as IRAS 17449+2320 with a strong magnetic field exhibit a very low stellar rotation rate.

Magnetospheric ejections
In previous studies of disc models that include a compact magnetosphere (  = 0.1) and a weak magnetic field, magnetospheric ejections have been found (Zanni & Ferreira 2013;Čemeljić 2019;Čemeljić et al. 2023).This type of coronal ejections can substantially modify the stellar rotation depending on the strength of the magnetic field and can produce a spin-down of the star ( Čemeljić et al. 2023).Remarkably, we find magnetospheric ejections in both disc models with a compact and with an extended magnetosphere (that is,   = 0.1 and   = 1.0, see for instance Fig. 2).We think that in the case of FS CMa-type stars the observed low-stellar rotation is produced by such magnetospheric ejections and, therefore, definitely deserves a detailed investigation which is beyond the scope of this study and should be studied elsewhere.

The mass flow
The mass flow rate crossing the whole domain is the integral of the mass flux inside and outside of the disc.It is given by Figure 4 shows the mass flow rate for four models (BSROT0 2a, BSROT0 2b, BSROT0 2c and BSROT0 2d), calculated at three different radius of the disc.The upper panel corresponds to  = 3.0 0 .We find that for the four models the rate of mass flux towards the central star is always negative or zero with values between zero to -0.3.This behaviour was expected since in that part there is not formation of strong outflows (see figure 1 bottom panel).Also this means that there is no inflow of gas towards the central star in the mid-plane of the disc.The middle panel in Fig. 4 corresponds to the value given in equation ( 23) truncation radius  = 3.97 0 .We can see some positive values in the models BSROT0 2a, BSROT0 2b and BSROT0 2c.In the bottom panel in Fig. 4 we show the case when  = 10.00 .Here we can see that for the models BSROT0 2c and BSROT0 2d the mass flux rate becomes positive after a time of  = 60 0 and  = 70 0 , respectively, which means that the accretion disc stops at a radius greater than  = 10.00 .

Backflow and radial velocity
Under some conditions, a phenomenon known as "backflow" occurs in the accretion disc (see Mishra et al. 2023, and references therein).This phenomenon appears when a part of the disc flow moves in the opposite direction, away from the central object and usually appears in the mid-plane of the disc.In the figure 5 we show the radial velocity for the models where the resistivity was obtained using the description given in Eq. ( 21) and considering the case when the rate of stellar rotation is  ★ = 0.2.
In the plots we show the temporal evolution of the radial velocity component calculated at three different radii.The blue circles correspond to  = 5 0 , the red triangles to  = 10 0 and the green stars symbols correspond to  = 15 0 .In the four models we can observe that when the radial velocity is calculated in a radius corresponding to  = 5 0 the radial velocity is approximately zero over time.However, when the radius is taken as  = 10 0 or  = 15 0 we find negative values for the radial velocity component.That means an intermittent backflow appers as a function of the time, i.e., that backflow persist for certain time intervals and disapear in others.
On the other hand, for a stellar rotation rate of  ★ = 0.1, considering a compact and extended magnetosphere (BSROT0 1c and BSROT0 1d models) we find again an intermittent backflow when we analyze the radial component of the velocity at radial distances from the star of  = 10 0 and  = 15 0 , as can be seen in the figure 6.It is important to point out that although in most of the BSROT0 1 and BSROT0 2 models the magnetic Prandtl number  m , is close to (or exceeds) the critical value of the Prandtl number ( crit m ≈ 0.6; see Mishra et al. 2023), we find intermittent backflow at the mid-plane of the disc.Therefore, our results support the argument given in Mishra et al. (2023) that the backflow is a physical effect rather than a numerical effect.

Observed disc structure in FS CMas stars
There are several observational studies in which the presence of thick disc structures around FS CMa-type stars has been suggested.For example, the analysis of FS CMa interferometric data by Hofmann et al. (2022) shows that the outer dust ring has an inclination angle of approximately 40 • .However, in the UV spectra a strong absorption of elements from the iron group is detected, indicating an enormous amount of gas at these angles.Among other common properties of these objects, the presence of expanding layers, which can slow down, has also been suggested (Kučerová et al. 2013).Clumps moving away from the star were detected in Jeřábková et al. (2016) and evidence of material infall in Kučerová et al. (2013).Another important effect was reported by Ellerbroek et al. (2015).Based in the interferometric observations of HD 50138 in  with AMBER.They found that the star HD 50138 is surrounded by an au-scale rotating gas disc and spherical halo.Among some of the results they reported are that the radius of the outer disc is 0.6au, the radius of the halo is 3au, and Δ in the halo is 60km/s.The results shown in our figure (7) can explain observed phenomena reported by Ellerbroek et al. (2015).The models BSROT0 1a and MSROT0 1a show that the structure of the disc is thicker and also we can observe a structure on the edge that could be explained as a halo surrounding the disc.
Another study that supports the idea that FS CMa-type stars are surrounded by thick disc was done by Zickgraf & Schulte-Ladbeck (1989).They found that the FS CMa stars are surrounded by a very extended and optically thick disc.Based on polarimetric observations at Calar-Alto, they found that the disc opening angle is ≈ 30 • for the stars MWC 645 and MWC 939.We find that in all our models there is an increase in the aspect ratio as it evolves over time.For instance, in the case of the BSROT0 1b model the aspect ratio increased 3.5 times for a corresponding time of  ★ = 130 stellar periods that means that in this case we found that the disc opening angle is ≈ 25 • , while for the BSROT0 2b model the disc opening angle is ≈ 33 • i.e. the aspect ratio increased 4.6 times for the same time period.In Fig. 1 one can see how the disc is thickened.

Magnetic field configuration inside the accretion disc
In the case of a compact magnetosphere (models with  m = 0.1), we find that the magnetic field configuration is far from the standard configuration of a vertical magnetic field (resulting from an initial dipole configuration).We find that the poloidal magnetic field lines are twisted (see for example Fig. 2).In such configuration, the magnetic reconection heating the corona is highly probable, that could explain the presence of the weak Raman lines found observationally in several FS CMa stars.Additionally, we can see in the upper panel in Fig. 2 that the magnetic field lines within the disc inflate and twist rapidly (within the first 10 orbital periods of the star), which gives rise to the formation of knots at different radial distances near the disc mid-plane.These knots substantially modify the density of the disc producing regions of low density.In other words, they can generate the formation of shallow gaps in the disc which can be displaced radially depending on the inflation rate of the magnetic field lines (see upper panel in Fig. 7).In addition, the knots can also move towards the corona region, dragging gas with them from the surface of the disc, so the disc suffers a thickening.It should be emphasized that regardless of the analytical form of the magnetic diffusivity described by equations (19-21), when  m = 0.1 we find the formation of at least a pair of knots tied to the mid-plane of the disc.This result implies the generation of strong turbulence within the disc which can modify the velocity components and the density distribution.
On the other hand, when  m = 1 we find that the formation of knots occurs in the corona region near the disc layer (see lower panel in Fig. 2), and therefore the perturbations in density and velocities within the disc are smaller.However, these perturbations produce a higher accretion rate in the mid-plane of the disc (see below) unlike other studies where strong magnetic fields are also considered (Miller & Stone 1997;Romanova et al. 2002;Čemeljić et al. 2023).We think that the main reason for obtaining a greater perturbation in density within the disc is due to the vertical dependence of the density.Since we use a vertical profile of the density that depends on the height through the polar -coordinate and, therefore the density of the disc decays for values of  close to the aspect ratio of the disc ℎ and facilities the inflation of the magnetic lines.Otherwise, as shown in Čemeljić et al. (2023) when a density-taper function ( =    −3/2 ) in the corona-disc interface region is used, the density becomes more noisy there without drastic density changes near the mid-plane since the magnetic field decreases by at least two orders of magnitude.

SUMMARY AND CONCLUSIONS
We have performed a set of 2.5 magnetohydrodynamical simulations to analyze the star-disc interaction of FS CMa type stellar objects that exhibit a strong magnetic field at the star surface ( ★ = 6.2;Korčáková et al. 2022).Based on the results of Paper I we set the value of the disc density  0 = 10 −13 gcm −3 , and we have considered as free parameters; the viscosity   , and the resistivity  m inside of the accretion disc, as well as the stellar rotation rate  ★ .By varying the resistivity, we analyze the effect of the magnetic field on the density of the disc and the corona region in the cases of a compact ( m = 0.1) and extended magnetosphere ( m = 1), respectively.
We find that a higher stellar rotation rate ( ★ = 0.2) can modify the density structure in the internal region of the disc, generating different layers when  ≤ 10 0 , which adds to the effect generated by a higher viscosity (  = 1) that produces the thickening of the accretion disc.In the case of a low-viscosity (  = 0.1), we find that in both models star-disc interaction with a compact or extended magnetosphere the fall of gas towards the star can be in the mid-plane.However, we do not rule out that the funnel effect may occur where the accretion of gas is conducted near the star's poles, since until the end of our simulations, the density structure in most of our models resembles a form similar to the funnel effect reported in less massive stars with a lower magnetic field (see for instance Romanova et al. 2002;Bessolaz et al. 2008;Zanni & Ferreira 2009, 2013;Čemeljić 2019).
On the other hand, since we consider a value of   = 0.1 < 0.685, we have analyzed if there is a backflow in the mid-plane of the disc which was previously reported (see Čemeljić 2019, and references therein).We find that regardless of the type of magnetosphere modeled there is intermittent backflow, as well as several knots forming in the magnetic field lines near and in the mid-plane of the disc.These knots produce perturbations in the density and velocity components, as well as the formation of shallow gaps in the density whose position depends on the inflation of the magnetic field lines.
Furthermore, we conclude that our findings may have direct implications for understanding the observational results for the study of post-mergers among FS CMa stars, because as we can see in our simulations the disc structure changes and becomes thicker significantly.Our simulations could also explain some other properties that have been found observationally in FS CMA stars (Kučerová et al. 2013), such as: material outflow and inflow, the formation of layers, and also the presence of the weak Raman lines found in several FS CMA stars than can be explained by corona heating due to the magnetic reconection.

Figure 1 .
Figure1.Upper panel shows logarithmic gas density for the BSROT 1b model with a low rate of stellar rotation  ★ ≡ Ω ★ /Ω br = 0.1.It can be seen the formation of the funnel effect as well as the formation of disc winds and magnetospheric ejections.Bottom panel shows logarithmic gas density for the BSROT 2b model with a higher rate of stellar rotation  ★ = 0.2.In this case we find the formation of multiple disc layers between  = 5 0 and  = 10 0 as well as the inflation of the disc that resembles the formation of a huge funnel.The white lines show the poloidal magnetic field lines.Note that in both models except for the increase in the stellar rotation rate in a factor 2, all other parameters remained fixed.

Figure 2 .
Figure 2. Upper panel shows logarithmic gas density for the BSROT 1a model with a low rate of stellar rotation  ★ = 0.1 as well as a low magnetic diffusivity   = 0.1.Bottom panel shows logarithmic gas density for the BSROT 1b model with a higher magnetic diffusivity   = 1.0.The white lines show the poloidal magnetic field lines.

Figure 3 .
Figure 3. Radial profiles at the midplane of the Sonic Mach number, angular velocity in the disc Ω (normalized by Ω ★ ) and the Keplerian velocity Ω K for the the BSROT 1a model calculated at  = 10 ★ .The green vertical line shows the position of the truncation radius given in Eq. (23).

Figure 5 .
Figure 5. Radial velocity in the mid-plane as a function of time, calculated at three different points  = 5 0 ,  = 10 0 and  = 15 0 .In all graphs presented has been considered that the rate of stellar rotation is  ★ = 0.2.The corresponding model is indicated at the top of each subplot.

Figure 6 .
Figure 6.Radial velocity at the mid-plane as a function of time, calculated at three different points  = 5 0 ,  = 10 0 and  = 15 0 .In all graphs presented has been considered that the rate of stellar rotation is  ★ = 0.1.At the top of each subplot the model to which it corresponds is indicated.
First column.Name of each model, the first letter indicates what type of magnetic diffusivity was used (B: Bessolaz, C: Classic and M: Modified, respectively).Second column.The stellar rotation rate.Third column.Rotation period of the star.Fourth column.The anomalous viscosity coefficient.Fifth column.Dimensionless coefficient of magnetic diffusivity.Sixth column.Prandtl number defined as  m ≡ 2  /3 m .

Figure 7 .
Figure 7.Comparison between different models.In this figure we show the results using a time period of 10 ★ , with a low stellar rotation rate  ★ = 0.1, and a low magnetic diffusivity   = 0.1 using the three different descriptions to obtain the resistivity discussed in section 2.5.The logarithmic gas density (color-scale background) changes from maximum value  = −0.1 in the disc to the minimum value  = −7.3 in the corona.The white lines show the poloidal magnetic field lines.

Table 1 .
Parameters of our numerical models.