© The Authors 2024

1. Introduction

Common envelope evolution (CEE) is a phase in the evolution of binary star systems where one star undergoes significant expansion and engulfs its companion within a shared envelope. The drag experienced by the companion star initiates its rapid spiral-in through the envelope, which leads to the transfer of energy and angular momentum to the envelope gas (Paczynski et al. 1976). CEE leads to two potential outcomes: either the stars merge into a single object or the rapid inspiral slows down. The exact causes of the inspiral deceleration remain elusive, but factors such as reduced gas density or corotation of gas with the companion likely play a major role (Roepke & De Marco 2022). Simulations suggest that the post-dynamical inspiral phase associated with weak drag from the shared envelope can persist for at least hundreds of orbits (e.g., Ricker & Taam 2012; Passy et al. 2012; Ohlmann et al. 2016a; Ivanova & Nandez 2016; Gagnier & Pejcha 2023). It is expected that such a weak drag gradually brings the binary closer together and ejects the envelope, ultimately resulting in a post-CEE binary (e.g., Ivanova et al. 2013; Clayton et al. 2017; Glanz & Perets 2018).

Accurately modeling CEE presents significant challenges due to the complex interplay of hydrodynamics, magnetohydrodynamics, radiative processes, and stellar evolution. While significant progress has been made through all-encompassing end-to-end simulations, interpreting these simulation results presents significant challenges due to a high degree of complexity resulting from a combination of multiple physical processes. A complementary approach, focusing on specific phases of CEE or distinct physical processes, is essential to achieve a deeper understanding. For instance,“wind tunnel” simulations conducted by MacLeod & Ramirez-Ruiz (2015), MacLeod et al. (2017), Cruz-Osorio & Rezzolla (2020), and De et al. (2020) characterized the flow properties around objects embedded within common envelopes to determine drag and accretion coefficients relevant for the rapid inspiral.

In Gagnier & Pejcha (2023), we carried out 3D hydrodynamical numerical simulations focusing exclusively on the post-dynamical inspiral phase of CEE. To have control over the simulation and to achieve a longer time step, we excised the central region containing the binary cores and we emulated the preceding dynamical plunge-in by depositing angular momentum in the envelope. This approach enabled us to perform a comprehensive analysis of the envelope structure and evolution, angular momentum transport, the short-term variability of accretion, and the estimation of the orbital contraction timescale. We found that the orbital contraction timescale is long, ∼10³ to ≳10⁵ orbits of the central binary, but it can become shorter than the thermal timescale of many envelopes. Furthermore, provided that there is gas remaining around the binary, the orbital contraction timescale never reaches zero, because the gas outside of the binary orbit cannot be kept in corotation and it develops spiral waves, which back-react on the orbit. Furthermore, our study of the variability of accretion exhibited remarkable similarities with what is seen in circumbinary disk simulations.

So far, the vast majority of CEE simulations have been performed only in the limit of hydrodynamics. It is possible that additional physical processes may introduce further complexity and potentially impact the results. In particular, in the field of stellar physics, magnetic fields play a central role in governing the dynamics of various systems, including protoplanetary and circumbinary disks. Within these disks, magnetic fields have a significant impact on crucial processes such as angular momentum transport (e.g., Papaloizou & Lin 1995; Balbus & Hawley 1998; Sano et al. 2004; Ji et al. 2006; Pessah et al. 2007; Shi et al. 2012), acceleration of winds (e.g., Bai & Stone 2013; Lesur 2021; Wafflard-Fernandez & Lesur 2023), and the launching of jets (e.g., Ferreira 1997; Gold et al. 2014; Qian et al. 2018; Vourellis et al. 2019; Saiki & Machida 2020). This raises the fundamental question of whether the influence of magnetic fields extends to CEE, where the highly turbulent environment could be conducive to their amplification. There exists a possibility that magnetic fields could significantly reshape the envelope’s structure and dynamics and could impact the binary’s orbital evolution. In fact, jet activity has been argued as a fundamental feature of CEE and shaping of planetary nebulae (e.g., Soker & Livio 1994; Soker 2016; Shiber et al. 2019; Hillel et al. 2022).

Due to numerical challenges of 3D magnetohydrodynamic (MHD) simulations of CEE, there are only a select number of papers addressing the potential effects of magnetic fields in this context. The first 3D MHD simulation of CEE was performed by Ohlmann et al. (2016b) who found that the amplification of magnetic energy was insufficient for magnetic fields to become dynamically significant during the dynamical inspiral phase of CEE. Schneider et al. (2019) conducted 3D MHD simulations of the merger of two massive stars, resulting in strong magnetic fields and yielding a rejuvenated merged star appearing younger and bluer than stars of its age. The most recent 3D MHD simulation of CEE was performed by Ondratschek et al. (2022). Their findings indicate that magnetic fields may have a substantial influence on shaping the morphology of emerging planetary nebulae by launching jet-like outflows from the immediate vicinity of the two central cores. Similar jet-like polar outflows were also observed in the 2D MHD post-common envelope circumbinary disk simulations of García-Segura et al. (2021). This scarcity of research presents an opportunity for conducting more thorough investigations into the role of magnetic fields in CEE.

In this paper, we aim to investigate the amplification of an initially weak seed magnetic field as it interacts with spiral density waves and hydrodynamical turbulence arising in the envelope from the time-changing gravitational force of the central binary system. We examined the influence of the resulting Lorentz force feedback on the structure and dynamics of the envelope and studied the impact on the orbital evolution of the central binary during the post-dynamical phase of CEE. Building upon our previous work in Gagnier & Pejcha (2023), we conducted 3D magnetohydrodynamical numerical simulations dedicated to the post-dynamical inspiral phase of the CEE. To emulate the outcome of the preceding dynamical inspiral phase, we injected angular momentum into the envelope following the methodology described in Morris & Podsiadlowski (2006, 2007, 2009), Hirai et al. (2021), and Gagnier & Pejcha (2023). We excised a central sphere containing the binary and we applied inner boundary conditions controlling the presence or absence of accretion. With these simulations, we performed a detailed analysis of the energy transfer both within and between the kinetic and magnetic energy reservoirs. Subsequently, we used similar techniques and diagnostics as those used by Gagnier & Pejcha (2023) to assess the impact of magnetic fields on the binary separation evolution timescale, the short-term variability of mass accretion onto the binary, the formation of overdensities, and angular momentum transport within the shared envelope.

This work is structured as follows: in Sect. 2, we introduce our physical model and describe the numerical setup used in our common envelope simulations. In Sect. 3, we present the results of our simulations. In particular, we study how kinetic and magnetic energy reservoirs are interconnected, and what contributes to their evolution. We then determine the scales at which the energy reservoirs interact using energy transfer analysis. By comparing our results with those of Gagnier & Pejcha (2023), we assess the impact of magnetic fields on the binary separation evolution timescale, on the short-term variability of mass accretion onto the binary, the formation of overdensities, and on angular momentum transport within the shared envelope. In Sect. 4, we summarize and discuss our results.

2. Physical model and numerical setup

We construct our post-dynamical inspiral model in the inertial frame at rest with the center of mass of the binary. We do not follow the previous evolution of the inspiraling binary, instead, we emulate its outcome following a procedure similar to Morris & Podsiadlowski (2006, 2007, 2009), Hirai et al. (2021), where the envelope is spun-up until a satisfactory amount of total angular momentum is injected. We describe the details of our implementation of this methodology in Gagnier & Pejcha (2023). We set the gravitational constant G, the total binary mass M = M₁ + M₂, the primary’s initial radius R, and the angular velocity $\sqrt{GM/R^3}$ to unity. The orbital velocity is fixed to ${\mathrm{\Omega }}_{\mathrm{orb}}=\sqrt{GM/a_{\mathrm{b}}^3}$, where a_b = r₁ + r₂ is the fixed binary separation, M₁ is the mass of the primary’s core located at {r₁, θ₁, φ₁}, and M₂ is the mass of the secondary object (either a main-sequence star or a compact object) located at {r₂, θ₂, φ₂}. The orbital period is P_orb = 2π/Ω_orb. The two objects are not resolved and are considered as point masses. To simplify our model and to make connection with our previous work, we consider an equal mass binary, q ≡ M₂/M₁ = 1, on a fixed circular orbit. The mass of the envelope is M_env = 2 in our units. Our choice of initial parameters for the binary and envelope broadly represents results of ab initio simulations for a range of progenitor types, as we detailed in Gagnier & Pejcha (2023). Because we are most concerned with the angular momentum transport within the common envelope in the two extreme regimes of mass and angular momentum accretion onto the binary rather than the specific details of the individual cores, we excise a central region encompassing the binary with a radius r_in = 0.625 a_b = R/10. In the rest of this Section, we describe the equations used for solving the problem (Sect. 2.1) and the initial conditions (Sect. 2.2). We then present our polar averaging implementation (Sect. 2.3) and the mesh structure (Sect. 2.4).

2.1. Equations of magnetohydrodynamics

We use Athena++ (Stone et al. 2020) to solve the equations of inviscid and ideal magnetohydrodynamics

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\rho \mathit{u}& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial \rho \mathit{u}}{\partial t}+\mathbf{\nabla }·(\rho \mathit{u}\otimes \mathit{u}-\mathit{B}\otimes \mathit{B}+P^{\ast }\mathit{I})& =-\rho \mathbf{\nabla }\mathrm{\Phi },\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial E}{\partial t}+\mathbf{\nabla }·((E+P\mathit{I})\mathit{u}-\mathit{B}(\mathit{B}·\mathit{u}))& =-\rho \mathbf{\nabla }\mathrm{\Phi }·\mathit{u},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}-\mathbf{\nabla }\times (\mathit{u}\times \mathit{B})& =0,\hfill \end{array}$$(4)

where P^*I is a diagonal tensor with components P^* = P + B²/2 = P + P_B, E = e + ρu²/2 + B²/2, e is the internal energy density, P = (Γ − 1)e, Γ = 5/3 is the adiabatic index, and Φ(r) is the gravitational potential of the binary,

$$\begin{array}{c}\hfill \mathrm{\Phi }(\mathit{r})=-{\displaystyle \sum _{i=1}^2\frac{G{M}_i}{|\mathit{r}-{\mathit{r}}_i|}.}\end{array}$$(5)

Heaviside–Lorentz units are used for electromagnetic quantities so that the magnetic permeability μ₀ = 1. In Athena++ (and in almost all higher order Godunov codes), the internal energy density is inferred from the difference between the total energy density E and the kinetic and magnetic energy densities. In the vicinity of the inner boundary, the plasma β parameter, β = P/P_B, may become very small locally, implying that the internal energy density is locally much smaller than the magnetic and kinetic energy densities. As a result, the MHD solver can locally return unphysically small or even negative internal energy densities. We circumvent this issue by solving the internal energy density equation (e.g., Stone & Norman 1992; Bryan et al. 2014; Takasao et al. 2015)

$$\begin{array}{c}\hfill \frac{\partial \check{e}}{\partial t}+\mathbf{\nabla }·(\mathit{u}\check{e})+(\mathrm{\Gamma }-1)\check{e}\mathbf{\nabla }·\mathit{u}=0,\end{array}$$(6)

in addition to the total energy density equation Eq. (3). Equation (6) ensures internal energy density positivity. The selection criterion operates on a cell by cell basis as follows

$$\begin{array}{c}\hfill e=\{\begin{array}{c}E-\rho u^2/2-B^2/2,E-\rho u^2/2-B^2/2<(1-{ϵ}_e)e\hfill \\ \check{e},\mathrm{otherwise}.\hfill \end{array}\end{array}$$(7)

We choose ϵ_e = 0.05. The system is completed by the same boundary conditions as in Gagnier & Pejcha (2023) for hydrodynamic variables to model the two regimes of accretion. Either the inner boundary supports the weight of the primary’s envelope preventing the fluid to flow through it, or the inner boundary is open to angular momentum and mass flow by imposing zero radial gradient of ρ, P, u_θ, and angular momentum in ghost zones. The horizontal components of the magnetic field are set to zero in inner ghost zones, and are copied from the last radial cell in the domain in outer ghost zones. The component normal to the boundaries is calculated by imposing the divergence-free constraint.

2.2. Initial conditions

As in Gagnier & Pejcha (2023), we assume that the gas in the envelope is initially in hydrostatic equilibrium and that it can be described by a polytropic equation of state ignoring the gas self-gravity and considering purely radial initial profiles. The radial density profiles are

$$\begin{array}{cc}\hfill \frac{\rho (r)}{\rho (r_{\mathrm{in}})}& ={(1+C(\frac{B}{3}(\frac{1}{r^3}-\frac{1}{r_{\mathrm{in}}^3})-A(\frac{1}{r}-\frac{1}{r_{\mathrm{in}}})))}^n,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill P(r)& =K{\rho }^{\mathrm{\Gamma }},\hfill \end{array}$$(9)

where n ≡ 1/(Γ − 1) = 3/2 is the polytropic index and

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill A^{\prime }& =A(\frac{1}{R}-\frac{1}{r_{\mathrm{in}}}),\hfill & \hfill B^{\prime }& =\frac{B}{3}(\frac{1}{R^3}-\frac{1}{r_{\mathrm{in}}^3}),\hfill \\ \hfill C& =\frac{\kappa -1}{B^{\prime }-A^{\prime }},\hfill & \hfill K& =\frac{(1-\mathrm{\Gamma })(B^{\prime }-A^{\prime })}{\mathrm{\Gamma }(\kappa -1)\rho {(r_{\mathrm{in}})}^{1/n}},\hfill \\ \hfill {\kappa }^n& =\rho (R)/\rho (r_{\mathrm{in}}).\hfill \end{array}\end{array}$$(10)

Density at the inner boundary ρ(r_in) is calculated from the prescribed total mass of the envelope. The primary star is initially embedded in a low-density medium to which we apply our outer diode-type boundary conditions and in which the envelope will expand later on. To model this low-density medium, we consider an atmosphere in hydrostatic equilibrium with constant ambient sound speed c_s, amb assuming $P=\rho c_{\text{s},\text{amb}}^2/\mathrm{\Gamma }$, which yields ${\rho }_{\text{ext}}\propto \exp \left(\Gamma /c_{\text{s},\text{amb}}^2\left(A/r-B/(3r^3)\right)\right)$. More details on the hydrodynamical initial conditions can be found in Gagnier & Pejcha (2023).

As is commonly done in (GR)MHD simulations of circumbinary disks (e.g., Noble et al. 2012; Shi et al. 2012; Lopez Armengol et al. 2021), the magnetic field is initialized as a single poloidal loop within the main body of the envelope at the end of the spin-up phase, with a vector potential A = {0, 0, A_φ} where

$$\begin{array}{c}\hfill A_{\phi }=A_{0}max(\rho -{\rho }_{\mathrm{cut}},0),\end{array}$$(11)

and ρ_cut = 10⁻⁵ρ_max. Here, ρ_max is the maximum density in the envelope at the moment when magnetic field is introduced and A₀ is computed to achieve the prescribed initial volume-averaged β parameter

$$\begin{array}{c}\hfill {\beta }_{\mathrm{m},\mathrm{i}}=\frac{{\int }_{\rho >{\rho }_{\mathrm{cut}}}P\mathrm{d}V}{{\int }_{\rho >{\rho }_{\mathrm{cut}}}{P}_B\mathrm{d}V}.\end{array}$$(12)

Specifically,

$$\begin{array}{c}\hfill A_0={\left(\frac{{\int }_{\rho >{\rho }_{\mathrm{cut}}}2P\mathrm{d}V}{{\beta }_{\mathrm{m},\mathrm{i}}{\int }_{\rho >{\rho }_{\mathrm{cut}}}\left||\mathbf{\nabla }\times max((\rho -{\rho }_{\mathrm{cut}}),0){\mathit{e}}_{\phi }\right|{|}^2\mathrm{d}V}\right)}^{1/2}.\end{array}$$(13)

In this work, we consider β_m, i = 10³ in the magnetized region, which corresponds to a field initially too weak to affect the dynamics of the envelope. In Fig. 1, we present a zoomed-in snapshot of density cross section in the xz plane as well as the poloidal magnetic field lines at the end of spin-up. In Table 1, we summarize parameters and outcomes of our runs.

Fig. 1. Zoomed-in snapshot of density cross section and initial poloidal magnetic field lines in the xz plane after the end of spin-up.

2.3. Polar averaging

The clustering of cells in the azimuthal direction near the polar axis makes the time step given by the Courant–Friedrichs–Lewy condition extremely small. To mitigate this issue, we follow Gagnier & Pejcha (2023) and use a polar averaging technique based on the Ring Average technique (Lyon et al. 2004; Zhang et al. 2019). The r- and θ-directed face-centered magnetic fields are also averaged in the azimuthal direction near the polar axis and the φ-directed field is updated from the Faraday’s law of induction, accounting for the electric field perturbations resulting from the reconstructed r- and θ-directed face-centered magnetic fields (e.g., Zhang et al. 2019). This procedure ensures ∇ ⋅ B = 0 within each cell and each averaged chunk of cells.

At some point in our simulations, the flow ends up stagnating near the polar axis, which leads to the formation of magnetic loops around the axis. If left unattended, the magnetic field can grow to large magnitude. To prevent this, we apply a resistive electric field along the axis after the reconstruction of the face-centered magnetic fields (see also Lyon et al. 2004). The resistive field is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_r^{\mathrm{axis}}{\mathrm{\Delta }}_r=\frac{cos{\theta }^{\mathrm{axis}}}{N_{\phi }{\tau }_{\mathrm{damp}}}{\displaystyle \sum _{k=1}^{N_{\phi }}{\mathrm{\Phi }}_{\phi }^B(k),}\end{array}$$(14)

where

$$\begin{array}{c}\hfill {\theta }^{\mathrm{axis}}=\{\begin{array}{cc}0\hfill & \mathrm{in}\mathrm{the}\mathrm{northern}\mathrm{hemisphere},\hfill \\ \pi \hfill & \mathrm{in}\mathrm{the}\mathrm{southern}\mathrm{hemisphere}.\hfill \end{array}\end{array}$$(15)

Here, N_φ is the number of cells in the φ direction adjacent to the polar axis at a given radius, Δ_r is the radial extent of the cell along the polar axis, ${\mathrm{\Phi }}_{\phi }^B(k)$ is the azimuthal magnetic flux through the kth cell face, and τ_damp is an arbitrary damping timescale which we choose to be a few time steps, 𝒪(10⁻⁵P_orb). By construction, the application of such a resistive electric field leaves ∇ ⋅ B unchanged. While primarily serving magnetic field stabilization, the application of such a resistive electric field may also result in the suppression of jets.

2.4. Mesh structure and convergence

At the root level, we set 512  ×  256  ×  256 active cells in {r, θ, φ}. The radial domain extends from r = 0.1 to r = 10 with geometric grid spacing so that the aspect ratio of the cells is approximately constant over the entire range of radial scales. The grid spacing in θ and φ directions is uniform with 0 < θ < π and 0 < φ < 2π. We use polar boundary conditions in the θ direction allowing free-flow through the pole (Zhu & Stone 2018; Stone et al. 2020), and periodic boundary conditions in the φ direction. On top of the root level, we add one level of adaptive refinement with criteria based on the azimuthal average of the second derivative error norm χ of a function σ (e.g., Lohner et al. 1987; Mignone et al. 2012). The norm χ is defined as

$$\begin{array}{c}\hfill \chi =\sqrt{\frac{\sum _d{|{\mathrm{\Delta }}_{d,+1/2}\sigma -{\mathrm{\Delta }}_{d,-1/2}\sigma |}^2}{\sum _d{\left(|{\mathrm{\Delta }}_{d,+1/2}\sigma |+|{\mathrm{\Delta }}_{d,-1/2}|+ϵ{\sigma }_{d,\mathrm{ref}}\right)}^2}}\ge {\chi }_r.\end{array}$$(16)

Here, Δ_r, ±1/2 = ±(σ_i ± 1 − σ_i) and σ_r, = |σ_i+1|+2|σ_i|+σ_i−1|. The value of the threshold χ_r depends on the specifics of the problem and on the chosen refinement variable σ. Here, we take σ = |B|.

We ensure that magnetorotational instability (MRI), which could occur in the simulation, is well resolved. We compute the quality factors in the three directions (e.g., Noble et al. 2010; Hawley et al. 2013) defined as

$$\begin{array}{c}\hfill Q_i=\frac{{\lambda }_{\mathrm{MRI},\mathrm{i}}}{{\mathrm{\Delta }}_i}=\frac{2\pi |v_{A,i}|}{\mathrm{\Omega }{\mathrm{\Delta }}_i},\end{array}$$(17)

where ${\mathit{v}}_{A,i}=\sqrt{B_i^2/\rho }$ is the i-component of the Alfvén speed in code units, and Δ_i is the size of individual cells in the i-direction. Using Q_i, we measure the resolvability R_i, which is the fraction of cells with a quality factor larger than that associated with marginal resolvability Q_min = 8 in the direction i (e.g., Sorathia et al. 2012; Parkin & Bicknell 2013). In Fig. 2, we show evolution of resolvability, volume-averaged quality factor, and the product Q_θQ_φ, for simulation run B. A value exceeding 200–250 for the product Q_θQ_φ is commonly regarded as a reliable indicator of convergence in MRI simulations (e.g., Narayan et al. 2012; Porth et al. 2019; Dhang et al. 2023). We find that while initially small due to limited envelope extent within the numerical domain, these quantities increase as the envelope expands and reaches the outer boundary, and surpass the threshold for convergence at t ≃ 55P_orb. This behavior can be attributed to the zero magnetic field amplitude we initially set in the low-density ambient medium encompassing the envelope. We conclude that our simulations have achieved numerical convergence.

Fig. 2. Temporal evolution of resolvability (panel a), quality factor (panel b), and QθQφ (panel c) in simulation run B. As the envelope expands within the numerical domain, these quantities progressively reach values that exceed the threshold values commonly associated with numerical convergence (horizontal dashed lines).

3. Results

We use a total of 4 million CPU hours on the Karolina cluster at IT4Innovations to perform our simulations. The parameters used for our simulation runs are outlined in Table 1. In Figs. 3 and 4, we present a zoomed-in snapshot of density and magnetic field amplitude cross section in the xz plane at t = 140P_orb for simulation B. In the rest of this section, we investigate the approach to quasi-steady state (Sect. 3.1), kinetic and magnetic energy budget (Sect. 3.2), evolution of the binary orbit (Sect. 3.3), angular momentum transport (Sect. 3.4), value of the α disk parameter (Sect. 3.5), time variability of accretion (Sect. 3.6), and behavior of the lump (Sect. 3.7).

Fig. 3. Zoomed-in snapshot of density cross section in the xz plane at t = 140Porb for simulation B. The texture, computed using the publicly available Line Integral Convolution Knit python package (Wafflard-Fernandez & Robert 2023), indicates the meridional streamlines.

Fig. 4. Zoomed-in snapshot of magnetic field amplitude cross section at t = 140Porb for simulation B. The texture indicates meridional magnetic field lines.

3.1. Quasi-steady state

After an initial transient period, the envelope settles into a quasi-steady state (QSS) during which the global properties of the envelope’s dynamics change only slowly. During this phase, matter continues to slowly accrete onto the central binary (if accretion is allowed) and simultaneously escapes from the numerical domain by the large-scale flow. To asses the turbulent fluxes of angular momentum and their quasi-steadiness, we apply Reynolds decomposition to the density, the gravitational potential of the binary, and the magnetic field

$$\begin{array}{cc}\hfill \rho & =\overline{\rho }+{\rho }^{\prime },\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \mathrm{\Phi }& =\overline{\mathrm{\Phi }}+{\mathrm{\Phi }}^{\prime },\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill B_i& =\overline{B_i}+B_i^{\prime },\hfill \end{array}$$(20)

where overlines indicate azimuthal averages of quantity q,

$$\begin{array}{c}\hfill \overline{q}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }q\mathrm{d}\phi .}\end{array}$$(21)

We also apply Favre decomposition to the velocity field (Favre 1965, 1969)

$$\begin{array}{c}\hfill u_i=\stackrel{\sim }{u_i}+u_i^{″},\end{array}$$(22)

where $\stackrel{\sim }{q}=\overline{\rho q}/\overline{\rho }$ is the density-weighted Reynolds average also known as Favre average. Taking the cross product of r with the momentum Eq. (2), multiplying by e_z and applying Reynolds average to the result yields the Reynolds averaged angular momentum evolution equation,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\partial \overline{\rho }s\stackrel{\sim }{u_{\phi }}}{\partial t}=-\mathbf{\nabla }·\left(s({\mathcal{W}}_{r\phi }+{\mathcal{W}}_{\theta \phi })\right)-\overline{\rho s\frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial \phi }},\end{array}\end{array}$$(23)

where s = r sin θ is the radial cylindrical coordinate,

$$\begin{array}{c}\hfill {\mathcal{W}}_{i\phi }=W_{i\phi }+w_{i\phi }\\ \hfill \end{array}$$(24)

is the total stress and

$$\begin{array}{cc}\hfill W_{i\phi }& =\overline{\rho }\stackrel{\sim }{u_i}\stackrel{\sim }{u_{\phi }}-\overline{B_i}\overline{B_{\phi }},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill w_{i\phi }& =\overline{\rho }\stackrel{\sim }{u_i^{″}{u}_{\phi }^{″}}-\overline{B_i^{\prime }{B}_{\phi }^{\prime }}.\hfill \end{array}$$(26)

We assume that a quasi-steady state is reached when the volume-averaged turbulent stress normalized to gas pressure (Shakura & Sunyaev 1973),

$$\begin{array}{c}\hfill {\alpha }_P=\frac{⟨w_{r\phi }⟩}{⟨P⟩},\end{array}$$(27)

and the normalized rφ–component of the Reynolds and Maxwell stresses, respectively

$$\begin{array}{cc}\hfill {\alpha }_K& =\frac{⟨2\rho u_r^{″}{u}_{\phi }^{″}⟩}{⟨\rho |{\mathit{u}}^{\mathbf{″}}{|}^2⟩},\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill {\alpha }_M& =-\frac{⟨2B_r^{\prime }{B}_{\phi }^{\prime }⟩}{⟨|{\mathit{B}}^{\mathbf{\prime }}{|}^2⟩},\hfill \end{array}$$(29)

become statistically time-independent (e.g., Simon et al. 2012; Parkin & Bicknell 2013). We show the time evolution of α_P, α_K and α_M in Fig. 5 for runs B and B′. We find that these three quantities reach quasi-steady values from t ≃ 100P_orb for simulation run B. α_K and α_M also reach quasi-steady values from t ≃ 100P_orb for simulation run B’, however, as α_K is small, it can change sign and lead to a substantial relative variation in α_P. Still, we assume that the system has reached a quasi-steady state from t ≃ 100P_orb for the two simulation runs. The time-averaged values of α_K, α_P and α_M at QSS are summarized in Table 1. The values of α_P and α_M at QSS for run B are in agreement with global and local magnetohydrodynamic simulations of accretion disks (e.g., Simon et al. 2012; Hawley et al. 2013; Parkin & Bicknell 2013) and indicate globally outward turbulent transport of angular momentum. These results, however, do not imply that the angular momentum turbulent flux cannot be locally directed inward nor that Maxwell stress dominates Reynolds stress. The radial and latitudinal dependence of angular momentum radial transport is investigated in Sect. 3.4.

Fig. 5. Evolution of the normalized volume-averaged stresses for simulation runs B (full lines) and B’ (dashed lines). We consider a quasi-steady state to be reached after t ≃ 100Porb.

3.2. Energy budget and dynamical relevance of magnetic fields

In this section, we determine how kinetic and magnetic energy reservoirs are interconnected, what contributes to their evolution, and the dynamical relevance of magnetic fields during the post-dynamical phase of CEE. In particular, we identify the source of magnetic energy amplification in Sects. 3.2.1 and 3.2.2, and determine the scales on which energy is transferred within and between energy reservoirs in Sect. 3.2.3. For the sake of brevity, most of this analysis is restricted to simulation run B, but results for B’ are analogous.

3.2.1. Magnetic energy budget

We show the evolution of the magnetic energy

$$\begin{array}{c}\hfill ⟨E_B⟩={\displaystyle \int }\frac{1}{2}\mathit{B}·\mathit{B}\mathrm{d}V\end{array}$$(30)

and kinetic energy

$$\begin{array}{c}\hfill ⟨E_K⟩={\displaystyle \int }\frac{1}{2}\rho \mathit{u}·\mathit{u}\mathrm{d}V,\end{array}$$(31)

for simulation B in Fig. 6. We find that, after a rapid amplification phase, the total magnetic energy saturates and then slowly decays. The initial rapid amplification of the magnetic energy results from the stretching, folding, and winding of the initial weak poloidal field by differential rotation and turbulence.

Fig. 6. Evolution of the total kinetic and magnetic energy and of their radial, latitudinal and azimuthal components for Gagnier & Pejcha (2023)’s simulation A (dashed lines) and this work’s simulation B (solid lines).

As the magnetic field strengthens, the Lorentz force begins to act back on the fluid, promoting more ordered and stable flow and leading to magnetic energy saturation. The magnetic tension associated with the strong toroidal field resulting from the winding up of the initial poloidal field, maintains the spiral density waves’ structure by resisting their deformation against radial perturbation. This can be seen in Fig. 7 where we compare density snapshots taken at t = 140P_orb for models A and B, revealing that model B exhibits sharper spiral waves and a more uniform density distribution compared to model A. Additionally, magnetic tension tends to align fluid motion with magnetic field lines, resulting in a preference for azimuthal kinetic energy over radial kinetic energy during the saturated phase, in contrast to nonmagnetic simulations (Fig. 6). Hence, the presence of magnetic fields, even if they are relatively weak, has a significant impact on the envelope’s structure and dynamics. This is discussed in more details later in this subsection, as well as in the following Sects. 3.2.2 and 3.2.3. Finally, the magnetic energy at saturation, approximately 10⁴⁴ erg, is similar to the findings from previous ab initio MHD simulations by Ohlmann et al. (2016b) and Ondratschek et al. (2022). However, in our case, the ratio between kinetic and magnetic energy is 10, which is significantly smaller than in their results, where it was roughly 1000.

Fig. 7. Zoomed-in snapshot of density cross section in the xy plane at t = 140Porb for Gagnier & Pejcha (2023)’s simulation A (top) and this work’s simulation B (bottom).

In Fig. 8, we show the kinetic and magnetic horizontally integrated energy spectra averaged over twenty spectra spanning three orbital periods between t = 137P_orb and t = 140P_orb during the saturated phase (see Appendix A and e.g., Baddour 2010; Parkin & Bicknell 2013). We find that kinetic energy is dominant on all scales with a diminishing ratio towards smallest scales where equipartition is approached. To assess the relevant physics associated with the amplification and saturation of the magnetic energy, we write the total magnetic energy evolution equation (see Appendix B),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ⟨{\dot{E}}_B⟩& ={\displaystyle \int (\mathit{B}\otimes \mathit{B}):\mathbf{\nabla }\mathit{u}\mathrm{d}V-\int E_B\mathbf{\nabla }·\mathit{u}\mathrm{d}V-\oint E_B\mathit{u}·\mathit{n}\mathrm{d}S,}\hfill \\ \hfill & =⟨{\dot{E}}_{B,\mathrm{stretch}}⟩+⟨{\dot{E}}_{B,\exp }⟩+⟨{\dot{E}}_{B,\mathrm{adv}}⟩,\hfill \end{array}\end{array}$$(32)

Fig. 8. Kinetic and magnetic horizontally integrated energy spectra during the saturation phase obtained with ℓmax = 128 for simulation run B. Opaque lines represent the spectra excluding ℓ = 0 components, while transparent lines include them. The oscillations visible at high k/kR when including ℓ = 0 components result from the our direct integration of Eq. (A.5) which involves spherical Bessel functions jℓ(kr) that are very oscillatory for kr ≫ ℓ.

where n is the outward-pointing unit vector at the domain boundaries surface, ⟨Ė_B,stretch⟩ is associated with the change in magnetic energy within the domain resulting from the stretching of the magnetic field lines by the fluid elements, ⟨Ė_B,exp⟩ corresponds to changes in magnetic energy due to expansion or compression effects, and ⟨Ė_B,adv⟩ corresponds to changes in magnetic energy from field advection through the domain boundaries.

The difference between the direct measurement of ⟨Ė_B⟩ and of the sum ⟨Ė_B,stretch⟩ + ⟨Ė_B,exp⟩ + ⟨Ė_B,adv⟩ corresponds to the numerical dissipation rate of magnetic energy by Joule heating ⟨Ė_B,diss⟩. Assuming numerical dissipation to be of the Ohmic form (e.g., Parkin & Bicknell 2013), we estimate the numerical resistivity as

$$\begin{array}{c}\hfill {\eta }_{\mathrm{num}}=\frac{⟨{\dot{E}}_{B,\mathrm{diss}}⟩}{⟨\mathit{B}\mathrm{\Delta }\mathit{B}⟩}=-\frac{⟨{\dot{E}}_{B,\mathrm{diss}}⟩}{\oint (\mathbf{\nabla }·\mathit{\sigma })·\mathit{n}\mathrm{d}S+\int {|\mathbf{\nabla }\times \mathit{B}|}^2\mathrm{d}V}.\end{array}$$(33)

Here, σ = B ⊗ B − E_BI is the total Maxwell stress tensor and we used integration by parts and the identity B × (∇×B) = ∇(B ⋅ B)/2 − (B ⋅ ∇)B. We show the different contributions to the evolution of the total magnetic energy density in Fig. 9 and our measure of the numerical resistivity η_num in Fig. 10. We find a magnetic Reynolds number ${\text{Re}}_M^{\text{QSS}}=\mathcal{O}({10}^3)$ at quasi-steady state for our models (see Table 1). However, we emphasize that the Ohmic dissipation formalism might not be valid and the measured resistivity likely depends on grid resolution, time integration method, spatial reconstruction scheme and Riemann solver (e.g., Rembiasz et al. 2017). In Fig. 9, we see that the dominant source of magnetic energy during both the growth and saturation phases is the stretching of the magnetic field lines by the velocity shear, which accounts for 95% of magnetic energy generation during QSS. The remaining 5% is produced by compression against magnetic pressure. The stretching of the magnetic field lines and the convergent transport of magnetic energy transfer kinetic energy into magnetic energy and effectively maintains dynamo action against dissipation and advection through the boundaries.

Fig. 9. Evolution of the rate of change of magnetic energy ⟨ĖB⟩ and of its components defined in Eq. (32), for simulation run B. The numerical dissipation rate of magnetic energy by Joule heating is measured as ⟨ĖB,diss⟩ = ⟨ĖB⟩ − ⟨ĖB,stretch⟩ − ⟨ĖB,exp⟩ − ⟨ĖB,adv⟩.

Fig. 10. Time evolution of the measured numerical resistivity (Eq. (33)) and viscosity (Eq. (40)) for simulation run B. At quasi-steady state, we find νnum ≃ 4.35 × 1015 cm2 s−1 and ηnum ≃ 7.11 × 1015 cm2 s−1.

3.2.2. Kinetic energy budget

To assess the relevant physics associated with the evolution of the kinetic energy budget and the dynamical relevance of magnetic fields in our simulations, we decompose the kinetic energy evolution equation into mean and turbulent contributions in Appendix C (see also Favre 1965, 1969). The Reynolds averaged mean kinetic energy evolution equation reads

$$\begin{array}{cc}\hfill \overline{{\dot{E}}_K^{\mathrm{mean}}(r,\theta )}& =\frac{1}{2}\frac{\partial \overline{\rho }(\stackrel{\mathbf{\sim }}{\mathit{u}}·\stackrel{\mathbf{\sim }}{\mathit{u}})}{\partial t}=-\mathbf{\nabla }·\left(\overline{\rho }\stackrel{\mathbf{\sim }}{\mathit{u}}\frac{(\stackrel{\mathbf{\sim }}{\mathit{u}}·\stackrel{\mathbf{\sim }}{\mathit{u}})}{2}\right)-\mathbf{\nabla }·(\overline{\rho }\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}\otimes {\mathit{u}}^{\mathbf{″}}}·\stackrel{\mathbf{\sim }}{\mathit{u}})\hfill \\ \hfill & +\overline{\rho }\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}\otimes {\mathit{u}}^{\mathbf{″}}}:\mathbf{\nabla }\stackrel{\mathbf{\sim }}{\mathit{u}}-\stackrel{\mathbf{\sim }}{\mathit{u}}·\mathbf{\nabla }\overline{P}+\stackrel{\mathbf{\sim }}{\mathit{u}}·(\mathbf{\nabla }·\overline{\mathit{\sigma }})-\stackrel{\mathbf{\sim }}{\mathit{u}}\overline{·\rho \mathbf{\nabla }\mathrm{\Phi }}+\stackrel{\mathbf{\sim }}{\mathit{u}}·(\mathbf{\nabla }·\overline{\mathit{\tau }}),\hfill \end{array}$$(34)

and the Reynolds averaged turbulent kinetic energy evolution equation reads

$$\begin{array}{c}\hfill \begin{array}{cc}& \overline{{\dot{E}}_K^{\mathrm{turb}}(r,\theta )}=\frac{1}{2}\frac{\partial \overline{\rho }(\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}·{\mathit{u}}^{\mathbf{″}}})}{\partial t}=-\mathbf{\nabla }·\left(\overline{\rho }\stackrel{\mathbf{\sim }}{\mathit{u}}\frac{(\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}·{\mathit{u}}^{\mathbf{″}}})}{2}\right)+\mathbf{\nabla }·(\overline{\mathit{\sigma }·{\mathit{u}}^{\mathbf{″}}}\hfill \\ \hfill & -\overline{\rho {\mathit{u}}^{\mathbf{″}}\frac{({\mathit{u}}^{\mathbf{″}}·{\mathit{u}}^{\mathbf{″}})}{2}}-\overline{P^{\prime }{\mathit{u}}^{\mathbf{″}}})-\overline{\rho }\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}\otimes {\mathit{u}}^{\mathbf{″}}}:\mathbf{\nabla }\stackrel{\mathbf{\sim }}{\mathit{u}}-\overline{\mathit{\sigma }:\mathbf{\nabla }{\mathit{u}}^{\mathbf{″}}}-\overline{{\mathit{u}}^{\mathbf{″}}}·\mathbf{\nabla }\overline{P}\hfill \\ \hfill & +\overline{P^{\prime }\mathbf{\nabla }·{\mathit{u}}^{\mathbf{″}}}-\overline{\rho {\mathit{u}}^{\mathbf{″}}·\mathbf{\nabla }{\mathrm{\Phi }}^{\prime }}+\overline{{\mathit{u}}^{\mathbf{″}}·(\mathbf{\nabla }·\mathit{\tau })},\hfill \end{array}\end{array}$$(35)

where τ = ν_numρc is the (numerical) viscous stress tensor, c = ∇u + (∇u)^⊺ − (2/3)(∇ ⋅ u)I is the shear tensor and ν_num is the numerical viscosity that we assume spatially constant. Finally, we integrate Eqs. (34) and (35) over the meridional plane, and multiply the result by 2π to obtain the mean and turbulent kinetic energy equations.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ⟨{\dot{E}}_K^{\mathrm{mean}}⟩& =⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,\mathcal{R}}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,\mathrm{turb}}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,P}^{\mathrm{mean}}⟩\hfill \\ \hfill & +⟨{\dot{E}}_{K,\mathcal{L}}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,\mathrm{\Phi }}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{mean}}⟩,\hfill \\ \hfill \end{array}\end{array}$$(36)

where

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{mean}}⟩& =-\langle \mathbf{\nabla }·\left(\overline{E_K^{\mathrm{mean}}}\stackrel{\mathbf{\sim }}{\mathit{u}}\right)\rangle ,\hfill \end{array}$$(37a)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathcal{R}}^{\mathrm{mean}}⟩& =-\langle \mathbf{\nabla }·(\overline{\rho }\stackrel{\mathbf{\sim }}{\mathit{u}}·\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}\otimes {\mathit{u}}^{\mathbf{″}}})\rangle ,\hfill \end{array}$$(37b)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{mean}}⟩& =\langle \overline{\rho }\stackrel{\sim }{{\mathit{u}}^{\mathbf{″}}\otimes {\mathit{u}}^{\mathbf{″}}}:\mathbf{\nabla }\stackrel{\mathbf{\sim }}{\mathit{u}}\rangle ,\hfill \end{array}$$(37c)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,P}^{\mathrm{mean}}⟩& =-\langle \stackrel{\mathbf{\sim }}{\mathit{u}}·\mathbf{\nabla }\overline{P}\rangle ,\hfill \end{array}$$(37d)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathcal{L}}^{\mathrm{mean}}⟩& =\langle \stackrel{\mathbf{\sim }}{\mathit{u}}·(\mathbf{\nabla }·\overline{\mathit{\sigma }})\rangle ,\hfill \end{array}$$(37e)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{\Phi }}^{\mathrm{mean}}⟩& =-\langle \stackrel{\mathbf{\sim }}{\mathit{u}}·\overline{\rho \mathbf{\nabla }\mathrm{\Phi }}\rangle ,\hfill \end{array}$$(37f)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{mean}}⟩& ={\nu }_{\mathrm{num}}\langle \stackrel{\mathbf{\sim }}{\mathit{u}}·(\mathbf{\nabla }·\overline{\rho \mathit{c}})\rangle .\hfill \end{array}$$(37g)

Here, $⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{mean}}⟩$ corresponds to the changes of mean kinetic energy from its advection through the domain boundaries, $⟨{\dot{E}}_{K,\mathcal{R}}^{\mathrm{mean}}⟩$ is associated with the transport of mean kinetic energy by Reynolds stress, $⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{mean}}⟩$ measures the destruction rate of mean kinetic energy into turbulence, $⟨{\dot{E}}_{K,P}^{\mathrm{mean}}⟩$ is the work done by pressure on the mean flow, $⟨{\dot{E}}_{K,\mathcal{L}}^{\mathrm{mean}}⟩$ is the work done by the Lorentz force on the mean flow, and $⟨{\dot{E}}_{K,\mathrm{\Phi }}^{\mathrm{mean}}⟩$ is the production rate of mean kinetic energy from the binary’s gravitational potential energy. The turbulent kinetic energy evolution equation in turn reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ⟨{\dot{E}}_K^{\mathrm{turb}}⟩& =⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{turb}}⟩+⟨{\dot{E}}_{K,\mathrm{trans}}^{\mathrm{turb}}⟩+⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{turb}}⟩+⟨{\dot{E}}_{K,P}^{\mathrm{turb}}⟩\hfill \\ \hfill & +⟨{\dot{E}}_{K,\sigma }^{\mathrm{turb}}⟩+⟨{\dot{E}}_{K,\exp }^{\mathrm{mean}}⟩,\hfill \end{array}\end{array}$$(38)

where

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{turb}}⟩& =-\langle \mathbf{\nabla }·\left(\overline{E_K^{\mathrm{turb}}}\stackrel{\mathbf{\sim }}{\mathit{u}}\right)\rangle ,\hfill \end{array}$$(39a)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{trans}}^{\mathrm{turb}}⟩& =\langle \mathbf{\nabla }·(\overline{\mathit{\sigma }·{\mathit{u}}^{\mathbf{″}}}-\overline{\rho {\mathit{u}}^{\mathbf{″}}\frac{({\mathit{u}}^{\mathbf{″}}·{\mathit{u}}^{\mathbf{″}})}{2}}-\overline{P^{\prime }{\mathit{u}}^{\mathbf{″}}})\rangle ,\hfill \end{array}$$(39b)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{turb}}⟩& =-⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{mean}}⟩,\hfill \end{array}$$(39c)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,P}^{\mathrm{turb}}⟩& =-\langle \overline{{\mathit{u}}^{\mathbf{″}}}·\mathbf{\nabla }\overline{P}\rangle ,\hfill \end{array}$$(39d)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\sigma }^{\mathrm{turb}}⟩& =-\langle \overline{\mathit{\sigma }:\mathbf{\nabla }{\mathit{u}}^{\mathbf{″}}}\rangle ,\hfill \end{array}$$(39e)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\exp }^{\mathrm{turb}}⟩& =\langle \overline{P^{\prime }\mathbf{\nabla }·{\mathit{u}}^{\mathbf{″}}}\rangle ,\hfill \end{array}$$(39f)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{\Phi }}^{\mathrm{turb}}⟩& =-\langle \overline{\rho {\mathit{u}}^{\mathbf{″}}·\mathbf{\nabla }{\mathrm{\Phi }}^{\prime }}\rangle ,\hfill \end{array}$$(39g)

$$\begin{array}{cc}\hfill ⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{turb}}⟩& ={\nu }_{\mathrm{num}}\langle \overline{{\mathit{u}}^{\mathbf{″}}·(\mathbf{\nabla }·\rho \mathit{c})}\rangle .\hfill \end{array}$$(39h)

Here, $⟨{\dot{E}}_{K,\mathrm{adv}}^{\mathrm{turb}}⟩$ corresponds to the changes in turbulent kinetic energy from its advection through the domain boundaries, $⟨{\dot{E}}_{K,\mathrm{trans}}^{\mathrm{turb}}⟩$ is the turbulent kinetic energy transport rate from the work done by the Lorentz and pressure forces and by turbulent velocity fluctuations, $⟨{\dot{E}}_{K,\mathrm{t}}^{\mathrm{turb}}⟩$ measures the production of turbulent kinetic energy from the interaction between mean shear and Reynolds stress, $⟨{\dot{E}}_{K,P}^{\mathrm{turb}}⟩$ is the work done by pressure on the turbulent flow, $⟨{\dot{E}}_{K,\sigma }^{\mathrm{turb}}⟩$ is the rate of production or destruction of turbulent kinetic energy from the interaction between turbulent shear and Maxwell stress, and $⟨{\dot{E}}_{K,\exp }^{\mathrm{mean}}⟩$ is the rate of change of turbulent kinetic energy due to pressure dilatation. Finally $⟨{\dot{E}}_{K,\mathrm{\Phi }}^{\mathrm{turb}}⟩$ is the production rate of turbulent kinetic energy from the work done by the multipole moments of the binary’s gravitational force on the turbulent flow. We note that because we apply Favre decomposition to the velocity field, the rate of change of the total kinetic energy is simply the sum of its mean and turbulent counterparts, that is, $⟨{\dot{E}}_K⟩=⟨{\dot{E}}_K^{\mathrm{mean}}⟩+⟨{\dot{E}}_K^{\mathrm{turb}}⟩$.

We show the evolution of the mean and turbulent kinetic energy change rates and of their components at QSS in Figs. 11 and 12. At QSS, we find that the binary orbital energy is the sole contributor to the envelope’s mean kinetic energy production, while other processes act as sinks of mean kinetic energy. The primary sinks for mean kinetic energy are the losses through the boundaries and the work done by pressure on the mean flow. Additionally, but to a lesser extent, mean kinetic energy is dissipated into turbulence and lost due to the advection of mean flow by the Reynolds stress and by the work done by the Lorentz force. Turbulent kinetic energy production is evenly split between the work done by the multipole moments of the binary’s gravitational force on the turbulent flow and the interaction between the mean shear and Reynolds stress. Conversely, the loss of turbulent kinetic energy mainly comes from its advection through the domain boundaries. Still, numerical dissipation, pressure dilatation, and interactions between turbulent shear and Maxwell stress significantly contribute to the loss of turbulent kinetic energy.

Fig. 11. Evolution of the rate of change of mean kinetic energy $⟨{\dot{E}}_K^{\mathrm{mean}}⟩$ and of its components defined in Eqs. (37a)–(37g), for simulation run B.

Fig. 12. Evolution of the rate of change of turbulent kinetic energy $⟨{\dot{E}}_K^{\mathrm{mean}}⟩$ and of its components defined in Eqs. (39a)–(39h), for simulation run B.

Similarly to the way we measured numerical Ohmic resistivity, we estimate the numerical kinematic viscosity, assuming it to be spatially constant, as

$$\begin{array}{c}\hfill {\nu }_{\mathrm{num}}=\frac{⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{mean}}⟩+⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{turb}}⟩}{\langle \stackrel{\mathbf{\sim }}{\mathit{u}}·(\mathbf{\nabla }·\overline{\rho \mathit{c}})\rangle +\langle \overline{{\mathit{u}}^{\mathbf{″}}·(\mathbf{\nabla }·\rho \mathit{c})}\rangle },\end{array}$$(40)

where $⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{mean}}⟩$ and $⟨{\dot{E}}_{K,\mathrm{diss}}^{\mathrm{turb}}⟩$ are obtained by substracting the various nondissipative components to the measured $⟨{\dot{E}}_K^{\mathrm{mean}}⟩$ and $⟨{\dot{E}}_K^{\mathrm{turb}}⟩$. We show our measurement of the numerical kinematic viscosity for simulation B in Fig. 10. We estimate the magnetic Prandtl number to be P_M ≃ 0.6, a value identical to that obtained by Parkin (2014) in the context of global simulations of magnetorotational turbulence using the PLUTO code (Mignone et al. 2012). However, astrophysical bodies are characterized by a wide range of magnetic Prandtl numbers which often substantially deviate from unity (e.g., Brandenburg & Subramanian 2005; Balbus & Henri 2008; Brandenburg 2009). Our measured P_M may thus not be realistic for common envelopes and cast doubts on the actual saturation values of α_M and α_P, that is, on angular momentum transport efficiency. Still, an asymptotic regime R_M > R_c in which α_P becomes P_M–independent may be reached in our simulations (e.g., Brandenburg 2009; Käpylä & Korpi 2011; Oishi & Mac Low 2011). Here, R_M = VL/η_num is the magnetic Reynolds number quantifying the relative importance of magnetic induction to magnetic diffusion, R_c is a critical Reynolds number of order 10³, and V and L are respectively a typical velocity and lengthscale of the flow. We reiterate here that our measure of the Prandtl number is only indicative and should be taken with caution. It is the result of strong approximations on the spatial uniformity of the numerical diffusivities and that the numerical diffusivity is equivalent to an Ohmic resistivity.

3.2.3. Energy transfer analysis

Here we analyze the scales at which the energy sources described in Sects. 3.2.1 and 3.2.2 operate and we focus attention to the scales at which magnetic energy is amplified. In order to do so, we perform a spectral energy transfer analysis, which is a commonly used method for understanding turbulent processes (e.g., Kraichnan 1967; Debliquy et al. 2005; Alexakis et al. 2007; Simon et al. 2009; Pietarila Graham et al. 2010; Lesur & Longaretti 2011; Rempel 2014; Grete et al. 2017, 2021; Du et al. 2023). We show the energy transfer equations and the resulting transfer functions in Appendix D. We emphasize that our approach involves horizontally integrated Fourier transformation (see Appendix A) for computing the transfer functions. Consequently, we are unable to assess the anisotropy of energy transfer in this analysis.

In Fig. 13, we show the energy transfer quantities averaged over 137 ≤ t/P_orb ≤ 140. In panel a, we display the magnetic energy transfer spectra associated with the stretching of magnetic field lines T_MT(k), the advection of magnetic energy within the magnetic energy reservoir T_MA(k), and with compression effects T_MC(k). In panel b, we show the spectra associated with the magnetic cascade T_MA(k)+0.5T_MC(k) and the kinetic to magnetic energy transfer T_MT(k)+0.5T_MC(k). Positive values of transfer functions at a given wavenumber k indicate gain of energy at a scale 2π/k. Conversely, negative values indicate a loss of energy at this scale. In panel a, we find that T_MT(k) is positive for all wavenumbers k, indicating that magnetic energy production occurs through line stretching across all scales. Conversely, T_MA(k) is only positive on scales smaller than ∼a_b/10, suggesting a forward magnetic energy cascade from larger to smaller scales through advection. Remarkably, both T_MT(k) and T_MA(k) peak at a scale of R_*/2 ≈ 3a_b, which is approximately the wavelength of the binary-driven spiral density waves, with almost equal but opposite values. This indicates that the bulk of magnetic energy production occurs at that scale through line stretching and that the newly generated magnetic energy is promptly advected towards the smallest scales. We find that the energy transfer rate by compression effects, T_MC(k), is overall weaker compared to T_MT(k) and T_MA(k). Unlike T_MT(k) and T_MA(k), transfer by compression effects T_MC(k) alternates between positive and negative values with varying k. This can be attributed to the compression–rarefaction pattern of spiral density waves. Although weaker, T_MC(k) plays a crucial role in the net magnetic energy production on the R_*/2 scale due to the near-perfect balance between T_MT(k) and T_MA(k) at this scale. In panel b, we see that, on scales larger than 2R_*, T_MT(k)+0.5T_MC(k) is in balance with T_MA(k)+0.5T_MC(k). Such balance indicates no net gain of magnetic energy at scales larger than 2R_* and thus implies the absence of large-scale magnetic field production. On smaller scales however, T_MT(k)+T_MA(k)+T_MC(k) is positive. Despite this, ⟨Ė_B⟩ ≲ 0 due to numerical dissipation of magnetic energy by Joule heating (see Figs. 6 and 9).

Fig. 13. Results of the energy transfer analysis. Panel a: Magnetic energy transfer as a function of the normalized wavenumber with kR = 2π/R* averaged over twenty spectra spanning three orbital periods during the saturated phase, 137 ≤ t/Porb ≤ 140. Positive values of transfer functions at a given wavenumber k (full line) indicate gain of energy at the scale 2π/k, negative values (dashed line) indicate a loss of energy at this scale. Panel b: Same as panel a, except splitting the magnetic energy transfer rate into magnetic cascade contribution TMA(k)+0.5TMC(k) and kinetic to magnetic energy transfer TMT(k)+0.5TMC(k).

In panel a of Fig. 13, we also show the time-averaged spectra associated with the gain and loss of kinetic energy from the work of the Lorentz force via magnetic tension and magnetic pressure, T_KL. The kinetic energy reservoir’s perspective offers a different view. Indeed, we see that kinetic energy is lost to magnetic energy on scales down to ∼0.4a_b. Approximately 20% of the energy drawn from the kinetic reservoir on these larger scales returns as kinetic energy on smaller scales. Remarkably, we find that most of the kinetic energy is lost to magnetic energy on the largest scales. However, the gain of magnetic energy on these scales is notably lower, which hinders the generation of significant large-scale magnetic fields during QSS. This highlights the nonlocality of interactions between magnetic and kinetic energy in spectral space (Alexakis et al. 2007). Unfortunately, we cannot pinpoint the specific origin or destination scales in spectral space. This would require a more complex shell-to-shell energy transfer analysis (e.g., Alexakis et al. 2007; Lesur & Longaretti 2011; Grete et al. 2017, 2021) that is beyond the scope of this work.

The lack of magnetic energy production on large radial scales implies that the amplified magnetic field may not be able to funnel and collimate a radial polar outflow, potentially impeding the formation of well-defined bipolar jets or jet-like outflows that often give rise to the characteristic bipolar shapes commonly observed in planetary nebulae (e.g., García-Segura et al. 1999, 2018, 2020; Zou et al. 2020; Ondratschek et al. 2022). Nonetheless, the presence of a strong toroidal magnetic field about the orbital plane introduces magnetic tension. This tension acts to decelerate the expansion of the envelope in directions perpendicular to the polar axis. Additionally, centrifugal forces, along with turbulent mixing, jointly contribute to the formation of intermittent, jittering low-density regions near the poles (see Fig. 3) which may facilitate the efficient channeling of outflows. Hence, the combined effects of magnetic tension, centrifugal forces, and turbulent mixing may ultimately contribute to the formation of nonspherical planetary nebulae.

3.3. Binary evolution and angular momentum conservation

Here, we address the impact of magnetic fields on the orbital evolution of the central binary. The commonly accepted model of common envelope evolution assumes a near-monotonic decrease in binary separation over time. In Gagnier & Pejcha (2023), we showed that this process is characterized by an orbital contraction timescale that can be much shorter than the thermal timescale of the envelope, particularly when the central binary is able to accrete mass and angular momentum from the surrounding envelope. In our setup, and as in Gagnier & Pejcha (2023), we fix the orbital parameters and predict the orbital evolution by measuring the rate of angular momentum transfer between the binary and the envelope without self-consistently considering any reciprocal influence from the envelope to the binary orbit. To assess the influence of magnetic fields on the secular evolution of binary separation, it is crucial to evaluate the torques acting within the shared envelope. These torques originate from both the quadrupolar moment of the gravitational potential and the angular momentum fluxes resulting from stresses across the domain boundaries. The conservation equation for angular momentum can be expressed as

$$\begin{array}{c}\hfill {\dot{J}}_z={\dot{J}}_{z,\mathrm{adv}}+{\dot{J}}_{z,\mathrm{grav}}+{\dot{J}}_{z,\mathrm{mag}},\end{array}$$(41)

where ${\dot{J}}_{z,\mathrm{adv}}$ is the advective torque associated with the loss of angular momentum through the boundaries,

$$\begin{array}{c}\hfill {\dot{J}}_{z,\mathrm{adv}}=-{\displaystyle {\oint }_{\partial R}\rho s{u}_{\phi }\mathit{u}·{\mathit{n}}_{\perp }\mathrm{d}S,}\end{array}$$(42)

${\dot{J}}_{z,\mathrm{grav}}$ is the gravitational torque exerted by the binary,

$$\begin{array}{c}\hfill {\dot{J}}_{z,\mathrm{grav}}=-{\displaystyle \int }\rho s\frac{\partial \mathrm{\Phi }}{\partial \phi }\mathrm{d}V,\end{array}$$(43)

and ${\dot{J}}_{z,\mathrm{mag}}$ is the magnetic torque,

$$\begin{array}{c}\hfill {\dot{J}}_{z,\mathrm{mag}}={\displaystyle {\oint }_{\partial R}s{B}_{\phi }\mathit{B}·{\mathit{n}}_{\perp }\mathrm{d}S.}\end{array}$$(44)

Here, n_⊥ is the outward-pointing unit vector at the boundaries’ surface. In Fig. 14, we show the evolution of these torques for runs B and B’. We find that at QSS, as in the nonmagnetic simulations of Gagnier & Pejcha (2023), the total angular momentum evolution is primarily driven by the outflow at the outer boundary, which is a consequence of envelope expansion and the finite radial extent of our numerical domain. Mass accretion hinders material buildup at the inner boundary and facilitates efficient horizontal turbulent mixing (as discussed in Gagnier & Pejcha 2023). Such efficient mixing weakens the gravitational torque, resulting in a modest transfer rate of angular momentum from the binary’s orbit to the envelope at QSS. Conversely, when accretion onto the binary is prevented (simulation run B’), material accumulation at the inner boundary and its stabilizing effect leads to an injection rate of angular momentum by the gravitational torque that is an order of magnitude larger than when accretion is allowed. The influence of magnetic fields on angular momentum evolution remains marginal in our two simulation runs. This is particularly evident in simulation run B’ where magnetic fields cannot be advected through the inner boundary by the fluid flow, unlike in simulation run B. Nevertheless, this impact very much depends on our selection of boundary conditions for the magnetic field.

Fig. 14. Moving average of the time evolution of the advective, magnetic, and gravitational torques for run B (panel a) and B’ (panel b). The measured total torque is indicated by a black line, positive values are denoted by solid lines and negative values by dashed lines.

As in Gagnier & Pejcha (2023), we set the orbital eccentricity e_b = 0 and enforce the binary mass ratio to q = 1. Consequently, we assume an equal distribution of mass and angular momentum accretion through the inner boundary between the two cores, resulting in $\dot{q}=0$. Additionally, we make the assumption that ė_b = 0. The time derivative of the binary’s angular momentum can then be expressed using the orbital separation evolution equation

$$\begin{array}{c}\hfill \frac{\dot{a_{\mathrm{b}}}}{a_{\mathrm{b}}}=\frac{\dot{M}}{M}(2\frac{M{\dot{J}}_{z,\mathrm{b}}}{\dot{M}{J}_{z,\mathrm{b}}}-3)=\frac{1}{{\tau }_{a_{\mathrm{b}}}},\end{array}$$(45)

where τ_ab is the orbital separation evolution timescale. If accretion onto the binary is disabled, Ṁ = 0 and ${\dot{J}}_{z,\mathrm{b}}=-{\dot{J}}_{z,\mathrm{grav}}$, otherwise $\dot{M}=-{\int }_{\partial R_{\mathrm{in}}}\rho u_r\mathrm{d}S$ and ${\dot{J}}_{z,b}=-{{\dot{J}}_{z,\mathrm{adv}}|}_{r=r_{\mathrm{in}}}-{\dot{J}}_{z,\mathrm{grav}}-{{\dot{J}}_{z,\mathrm{mag}}|}_{r=r_{\mathrm{in}}}$. We show the orbital separation evolution timescale for our two MHD runs B and B’, along with the non-MHD runs A and A’ from Gagnier & Pejcha (2023) in Fig. 15. As expected from the previous section, magnetic fields play little to no role in the binary separation evolution. Hence, as in the nonmagnetic simulations of Gagnier & Pejcha (2023), we find that the orbital separation evolution timescale reaches a quasi-steady value of the order 10–100 yr when accretion onto the binary is allowed, and increases with time as

$$\begin{array}{c}\hfill {\tau }_{a_{\mathrm{b}}}\propto {10}^{\gamma t/P_{\mathrm{orb}}},\end{array}$$(46)

Fig. 15. Time evolution of the orbital separation evolution timescale for our two MHD runs B and B’, along with the non-MHD runs A and A’ from Gagnier & Pejcha (2023). The black dashed line represents τab = 50 yr, while the black dashdotted line corresponds to τab = λ10γt/Porb with λ = 0.42 and γ = 0.0115.

with γ ≃ 0.0115, when accretion onto the binary is prevented.

3.4. Angular momentum transport

Radial transport of angular momentum within the common envelope can either facilitate or hinder the binary’s post-dynamical inspiral through the generation of accretion flows. This process also shapes the envelope’s structure by introducing density distribution asymmetries that enable angular momentum exchange between the binary and the envelope via the gravitational torque. These interactions further contribute to the shaping of emerging planetary nebulae resulting from the common envelope outflow. Magnetic fields are also well known to enhance or impede angular momentum transport within various astrophysical systems, as well as to generate torques, instabilities, and outflows, such as winds or jets. Here, we analyze angular momentum transport within the shared envelope, evaluating the morphology and strength of this process while assessing the influence of each physical mechanism.

In Fig. 16, we show the radial profile of the gravitational torque and mean-flow and turbulent contributions to the local angular momentum transfer rate across the common envelope, averaged over 140 ≤ t/P_orb ≤ 160 for simulations B (panel a) and B’ (panel b). The transfer rates are defined as (see Eq. (23))

$$\begin{array}{cc}\hfill {\overline{\dot{J}}}_{z,\mathrm{adv}}(r,t)& =-{\displaystyle {\oint }_{\partial r}s\overline{\rho }{\stackrel{\sim }{u}}_{\phi }{\stackrel{\sim }{u}}_r\mathrm{d}S,}\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill {\overline{\dot{J}}}_{z,\mathcal{R}}(r,t)& =-{\displaystyle {\oint }_{\partial r}s\overline{\rho }\stackrel{\sim }{u_{\phi }^{″}{u}_r^{″}}\mathrm{d}S,}\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill {\overline{\dot{J}}}_{z,\mathrm{mag}}(r,t)& ={\displaystyle {\oint }_{\partial r}s{\overline{B}}_{\phi }{\overline{B}}_r\mathrm{d}S,}\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill {\overline{\dot{J}}}_{z,\sigma }(r,t)& ={\displaystyle {\oint }_{\partial r}s\overline{B_{\phi }^{\prime }{B}_r^{\prime }}\mathrm{d}S,}\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill {\overline{\dot{J}}}_{z,\mathrm{grav}}(r,t)& ={\displaystyle {\int }_r^{R_{\mathrm{domain}}}\left({\int }_{\partial r}s\overline{\rho \frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial \phi }}\mathrm{d}S\right)\mathrm{d}r.}\hfill \end{array}$$(51)

Fig. 16. Gravitational torque and mean-flow and turbulent contributions to the local angular momentum transfer rate across the common envelope for model B (panel a) and B’ (panel b), averaged in time from t = 140Porb to t = 160Porb. Dashed lines depict negative values, whereas solid lines represent positive values.

When accretion onto the binary is allowed (simulation run B), we find that angular momentum is transported inward for r ≲ 1.4R_*, and outward further out in the envelope. It is dominated by the contribution from the mean-flow advection ${\overline{\dot{J}}}_{z,\mathrm{adv}}(r,t)$ which is directed inward for r ≲ 1.8R_*, and outward further out. The transport of angular momentum associated with the mean magnetic field ${\overline{\dot{J}}}_{z,\mathrm{mag}}(r,t)$ is directed outward across the envelope, with negligible amplitude. Within r ≲ 0.75R_*, turbulent Reynolds ${\overline{\dot{J}}}_{z,\mathcal{R}}(r,t)$ and Maxwell stresses ${\overline{\dot{J}}}_{z,\sigma }(r,t)$ contribute equally and significantly to outward angular momentum transport. Moving outward, the Reynolds stress changes direction intermittently, while the turbulent Maxwell stress maintains an outward direction with comparable magnitude. Overall, we find the combined effect of turbulent Reynolds and Maxwell stresses on radial transport of angular momentum to be significant. The combined effect damps the inward transport of angular momentum driven by the mean flow by a factor ∼2 for r ≲ 1.4R_*, it reverses the direction of angular momentum transport for 1.4 ≲ r/R_* ≲ 1.8, and the turbulent Maxwell stress slightly enhances outward transport for 1.8 ≲ r/R_* ≲ 4. Conversely, the binary’s gravitational torque impact on the radial angular momentum transport is negligible throughout the envelope. However, a notable feature is the presence of global radial oscillations about zero, which arise due to the gravitational coupling between the binary potential and the spiral density waves originating near the binary and propagating out in the envelope (Cimerman & Rafikov 2024).

When accretion onto the binary is prevented (simulation run B’), angular momentum is transported outward across the envelope. Notably, the positive values of ${\overline{\dot{J}}}_z$ within the very inner envelope result from the injection of angular momentum by the binary’s gravitational torque and they do not indicate inward transport. We further note that this injection of angular momentum is more pronounced in simulation B’ compared to simulation B owing to material accumulation at the reflecting inner boundary. In simulation B’, and contrary to simulation B, magnetic fields have a negligible impact on radial angular momentum transport. The primary driver remains advection by the mean flow, with the Reynolds stress also making a significant contribution. The Reynolds stress alternates between dampening and enhancing outward transport across varying radii from the central binary.

Because the time-periodic gravitational body force exerted by the binary system on its surrounding envelope is zero at the poles and maximum in the orbital plane, the dynamics of common envelopes is highly dependent on latitude. It is therefore essential to analyze the latitudinal dependence of angular momentum transport. To do so, we first show in Figs. 17 and 18 the total angular momentum radial advective flux F_adv = s𝒲_rφ (see Eq. (23)). We see that in simulations B and B’ and as in the nonmagnetic simulations of Gagnier & Pejcha (2023), the angular momentum advective flux is directed outward in an equatorial disk-like structure. Above and below such disk, the angular momentum advective flux is directed inward. We assess the relative significance and latitudinal variation of each component of the total angular momentum radial advective flux in Figs. 19 and 20. The top panel shows that the morphology of the angular momentum radial flux is dominated by the nonmagnetic component. The radial transport of angular momentum by magnetic fields is directed outward at all latitudes and radii for simulation B. Consequently, it enhances outward transport in the equatorial disk-like structure, and damps the inward transport elsewhere. The morphology of magnetic fields driving radial transport of angular momentum is slightly more complex in simulation B’ as its direction and amplitude vary with r and θ. From the second and third panels of Figs. 19 and 20 showing the mean and turbulent components of the total flux, we see that the turbulent transport associated with the Reynolds stress is less efficient than transport by the mean flow. This result is similar to the nonmagnetic simulations of Gagnier & Pejcha (2023). Still, the complex morphology and nonnegligible amplitude of the turbulent transport may significantly damp or enhance the angular momentum radial transport locally and may even reverse the latitudinally integrated angular momentum radial transport. We see that the magnetic radial transport of angular momentum is dominated by the turbulent Maxwell stress while the transport by mean (axisymmetric) field is negligible.

Fig. 17. Reynolds average of the total radial angular momentum flux s𝒲rφ, averaged over twenty orbital periods for 140 ≤ t/Porb ≤ 160 for simulation B.

Fig. 18. Same as Fig. 17 but for simulation B’.

Fig. 19. Reynolds average of the mean and turbulent contribution of the of the total radial angular momentum flux (see Eqs. (25) and (26)), averaged over twenty orbital periods for 140 ≤ t/Porb ≤ 160 for simulation B.

Fig. 20. Same as Fig. 19 but for simulation B’.

3.5. The value of the disk α parameter

Recently, Tuna & Metzger (2023) have developed an α-disk model for the long-term evolution of post-common envelope circumbinary disks. To check their assumptions and to facilitate similar studies in the future, it is important to constrain the value of α more accurately. In our analysis, presented in Sect. 3.1 and summarized in Table 1, we measured the volume-averaged turbulent stress normalized to the gas pressure, α_P, as well as the normalized rφ component of the Reynolds and Maxwell stresses. These measurements indicate that at QSS, angular momentum is globally transported outward when accretion is not prevented. The values we obtain are in agreement with previous global and local magnetohydrodynamic simulations of accretion disks (e.g., Simon et al. 2009; Hawley et al. 2013; Parkin & Bicknell 2013).

However, while many simulations of circumbinary disks are dominated by MRI with α_P mainly driven by the Maxwell stress that is predominantly positive throughout the disk (e.g., Shi et al. 2012; Shi & Krolik 2015), our simulations reveal a contrasting behavior. Magnetic fields in our accreting simulations are weak compared to the strength of the mean flow driven by the binary and have consequently limited direct dynamical impact (in agreement with Ondratschek et al. 2022). Yet, as we have seen in Sect. 3.2.1, magnetic fields have a considerable effect on the envelope’s density structure and therefore indirectly on angular momentum transport through the mean flow and through Reynolds stress. The Reynolds stress is the main contributor to the local turbulent flux of angular momentum. However, because it exhibits significant variations in sign and amplitude throughout the shared envelope (see also Gagnier & Pejcha 2023), its contribution to α_P is of the same order as that of the turbulent Maxwell stress. The turbulent Maxwell stress enhances the outward and damps the inward turbulent transport of angular momentum without changing its global morphology. We obtain similar results when accretion through the inner boundary is prevented by reflecting boundary conditions (Figs. 18 and 20). Angular momentum transport is dominated by the mean flow and consists of an outward transport in a disk-like structure and an inward transport along the polar axis. The main difference comes from the more complex morphology of the Maxwell and Reynolds stress and their overall weaker amplitude resulting in near-zero net radial turbulent transport of angular momentum. The global morphology of the turbulent transport of angular momentum is rather complex as its direction varies with latitude and distance to the central binary (see Figs. 21 and 22). Hence the relevance of employing a spatially constant value of α_P in post-CE circumbinary disks models may be questionable.

Fig. 21. Reynolds-averaged turbulent stress normalized to gas pressure, αP, averaged over twenty orbital periods for 140 ≤ t/Porb ≤ 160 for simulation B.

Fig. 22. Same as Fig. 21 but for simulation B’.

Our findings also reveal that the radial transport of angular momentum is predominantly driven by the mean flow, which essentially consists of spiral density waves. This poses a significant challenge to the use of the conventional α–viscosity formalism in 1D models, as it cannot directly account for the nondiffusive nature of angular momentum transport associated with these waves. Alternative models are needed to capture the intricate dynamics and angular momentum transport involving spiral density waves in post-CE circumbinary disks.

3.6. Time variability of accretion

In Fig. 23, we show a space-time diagram of the mass flux crossing the inner boundary normalized by its maximum value in the time interval 140 ≤ t/P_orb ≤ 160 for our accreting model B. As in the nonmagnetic simulations of Gagnier & Pejcha (2023), we see that mass accretion exhibits periodic variability at all colatitudes. On the orbital plane, a variability frequency of 2Ω_orb is clearly visible. This frequency is associated with the quadrupolar moment contribution to the binary potential. In addition to this expected variability frequency, Gagnier & Pejcha (2023) further found an additional low-frequency modulation of accretion ω ≃ 0.2Ω_orb, which they attribute to the presence of an eccentric and tilted overdensity, successively feeding the individual binary components through accretion streams. Remarkably, such mass accretion variability frequency are identical to those measured in circumbinary disks simulations where it is attributed to the orbiting frequency of a nonaxisymmetric overdensity (or “lump”). We show the power spectral density of the total mass accretion rate ṁ, as well as the mass accretion rates measured in various opening angles about the orbital plane in Fig. 23, panel c. We find the 0.2Ω_orb frequency is present, suggesting that magnetic fields do not prevent the formation of the overdensity observed in the nonmagnetic simulations of Gagnier & Pejcha (2023). It is however likely that higher magnetization would weaken the lump and its induced low-frequency accretion modulation (Noble et al. 2021). In addition to the 0.2Ω_orb frequency, we find a clear and dominating peak at ω ≃ 0.1Ω_orb. We investigate its origin in Sect. 3.7. We note that the 0.1Ω_orb and 0.2Ω_orb frequencies are not visible in the power spectral density of the total mass accretion rate because of the asynchronicity of mass accretion between colatitudes. These results suggest that local analyses are necessary when studying short-term variability of accretion in CEE (Gagnier & Pejcha 2023).

Fig. 23. Detailed view on the variability of mass flux for model B. Panel a: Space-time diagram of the local mass flux through the inner boundary. Panel b: Time evolution of the mass accretion rate onto the binary. Panel c: Power spectral density of the total mass accretion rate onto the binary (black line) and of the mass accretion rate in various colatitude ranges about the orbital plane (colored lines). Green lines correspond to the range π/4 ≤ θ ≤ 3π/4, blue lines to the range 3π/8 ≤ θ ≤ 5π/8, and orange lines to the range 7π/16 ≤ θ ≤ 9π/16.

3.7. The lump

In the context of circumbinary disks, the appearance of a lump is commonly thought to be the result of the interaction between gas streams inside the cavity and the cavity edge (e.g., Shi et al. 2012; Noble et al. 2012). However, recent work by Mignon-Risse et al. (2023) and Rabago et al. (2023) suggest an alternative lump formation channel consisting of the merging of Rossby Wave Instability (RWI; Lovelace et al. 1999; Li et al. 2000) induced large-scale vortices into a single lump. Cimerman & Rafikov (2024) further show that such large-scale vortices likely also play a major role in angular momentum transport in CBDs by launching vortex-driven spiral density waves (see also Huang et al. 2019). In circumbinary disks, RWI can be triggered from an axisymmetric bump in the inverse vortensity associated with the strong radial gradient of density at the cavity edge. In common envelopes however, there is likely no cavity, and such strong vortensity (or magnetized vortensity) maxima may not naturally emerge. Still, Rossby wave instability can be triggered by the vortensity production caused by the nonlinear damping of spiral density waves (Coleman et al. 2022). Here, we investigate the existence of nonaxisymmetric overdensities in our simulations and we explore their behavior and origin.

In Fig. 24, we show a snapshot of the relative surface density perturbation $(\mathrm{\Sigma }-\overline{\mathrm{\Sigma }})/\overline{\mathrm{\Sigma }}$ at t = 140P_orb for simulation run B. The surface density is defined as

$$\begin{array}{c}\hfill \mathrm{\Sigma }={\displaystyle {\int }_{7\pi /16}^{9\pi /16}\rho rsin\theta \mathrm{d}\theta .}\end{array}$$(52)

Fig. 24. Snapshot of relative surface density perturbation $(\mathrm{\Sigma }-\overline{\mathrm{\Sigma }})/\overline{\mathrm{\Sigma }}$ in the xy plane, at t = 140Porb.

We find strong and structured nonaxisymmetric overdensities taking the form of spiral density waves with a large pitch angle. The presence of a single well-defined spiral arm reaching the binary in Fig. 24 suggests that m = 1 modes also play a major role in the shaping of overdensities. In Fig. 25, we investigate the origin of such overdensities by first showing the nonaxisymmetric density perturbation on the orbital plane as well as the vortensity perturbation early in our simulation at t = 25P_orb. In the inner part of the envelope, we find small scale turbulence to be fully developed, which is likely a consequence of stratified shear instability emerging from the interaction between spiral density waves generated by the binary’s torque and the background velocity shear. We further note the presence of four large anticyclonic vortices located at r ≃ 1. These vortices likely emerge from Rossby Wave Instability, triggered by the vortensity production at this radius due to the nonlinear damping of spiral density waves (Coleman et al. 2022). The propagation of spiral density wavefronts through such vortices induces their deformation, as is evident in panel b. This deformation of the wavefront entails the emergence of crests where the deformed density waves converge. Notably, two out of the four vortices exhibit sufficient strength to cause the associated wavefront crests to give rise to the formation of persistent, slowly orbiting overdensities taking the form of m = 2 spiral density waves with large pitch angle. The existence of such overdensities prevent the spiral density waves generated by the binary’s torque from propagating freely outward, resulting in material accumulation and the formation of a nonaxisymmetric lump (see Gagnier & Pejcha 2023, for more details).

Fig. 25. Snapshot of density perturbation $(\rho -\overline{\rho })/\overline{\rho }$ (panel a) and of vortensity deviation $\zeta -\overline{\zeta }$, where ζ = ρ−1∇ × u, in the rφ plane at t = 25Porb.

Eventually, the inner parts of the two spiral arms become sufficiently dense for the gravitational force exerted on them by the central binary to overcome the pressure gradient driving the envelope’s expansion, leading to their inward migration. In Fig. 26, we see that the two overdensities exhibit slightly different inward migration speeds due to the imperfect mass distribution symmetry. Such inward migration of overdensities within the expanding envelope induces strong and unstable horizontal shear, giving rise to two prominent vortices. As the faster-moving overdensity approaches the binary, it experiences substantial torque and ultimately reaches the inner boundary as a singular accretion stream. This event locally disrupts the m = 2 symmetry, favoring instead m = 1 symmetry in the inner part of the envelope (Fig. 27). Over time, the density of this single spiral arm diminishes as material is accreted by the binary and is eventually destroyed by the m = 2 spiral density waves and by tidal disruption of the overdense stream, thereby restoring the m = 2 symmetry in the binary’s vicinity. This cycle repeats every ∼10P_orb and it therefore associated with the ω ≃ 0.1Ω_orb frequency observed in Fig. 23. We note that this process and thus this characteristic frequency were not present in Gagnier & Pejcha (2023)’s nonmagnetic runs. In simulation B, the formation and persistence of m = 1 accretion streams can be attributed to the stabilizing effect of the toroidal magnetic fields preventing the destruction of overdensities through magnetic tension, allowing these structures to progressively increase in density. As these overdensities become denser, their gravitational attraction to the inner binary system grows stronger, leading to their inward migration. Their relatively high density enables them to reach the inner boundary and feed the binary before being destroyed.

Fig. 26. Same as Fig. 25 but at t = 45Porb.

Fig. 27. Same as Figs. 25 and 26 but at t = 140Porb.

4. Discussions and conclusions

In this paper, we have significantly extended our previous work in Gagnier & Pejcha (2023) on the post-dynamical inspiral phase of CEE by including magnetic fields in 3D simulations. Our first aim was to determine the sources of magnetic energy amplification, and the relative importance of magnetic to kinetic energy reservoir size. We found that magnetic energy amplification arises primarily from the stretching, folding, and winding of the initial weak poloidal field due to differential rotation and turbulence. As the magnetic field strengthens, it stabilizes fluid motion, favoring azimuthal kinetic energy over radial kinetic energy during saturation. Magnetic fields significantly impact the envelope’s structure and dynamics, with the magnetic energy reaching levels similar to the previous MHD simulations of Ohlmann et al. (2016b) and Ondratschek et al. (2022), but with a much lower kinetic-to-magnetic energy ratio. Energy spectra show kinetic energy dominance at all scales, but decreasing towards smaller scales. Our analysis identifies stretching of magnetic field lines by velocity shear as the dominant source of magnetic energy, contributing 95% during quasi-steady state with the remaining 5% resulting from compression against magnetic pressure.

Our second aim was to determine how kinetic and magnetic energy reservoirs are interconnected and what contributes to their evolution during the post-dynamical phase of CEE. We first decomposed the kinetic energy evolution equation into mean and turbulent contributions to determine the relative contribution of the different processes. We found the impact of magnetic fields on the mean kinetic energy evolution to be negligible, and thus not to directly affect the mean flow. However, the interaction between turbulent shear and Maxwell stress leads to a nonnegligible sink of turbulent kinetic energy. We then employed energy transfer analysis to investigate the scales at which energy is transferred between kinetic and magnetic energy reservoirs. This analysis revealed that magnetic energy production by the stretching of magnetic field lines occurs at all scales and that the advection of magnetic energy within the magnetic energy reservoir indicates a forward cascade from larger to smaller scales. Notably, both field line stretching and magnetic advection peak at a critical scale of approximately 3a_b, that is, approximately the wavelength of the binary-driven spiral density waves. Compression effects, though weaker, play a significant role in net magnetic energy production due to the balance between stretching and advection at this scale. Our analysis also shows that the scales at which kinetic energy is lost to magnetic energy and the scales at which magnetic energy is received from kinetic energy are different. This highlights the nonlocality of interactions between magnetic and kinetic energy in spectral space.

Our analysis indicates the absence of magnetic energy production on large radial scales. This suggests potential challenges in funneling and shaping radial polar outflows in planetary nebulae, possibly hindering the formation of well-defined bipolar structures. However, the presence of a strong toroidal magnetic field on the orbital plane, coupled with centrifugal forces and turbulent mixing, could respectively slow down envelope equatorial expansion and channel jittery and irregular outflows near the poles, thus leading to the emergence of highly nonspherical planetary nebulae.

By comparing our results to the nonmagnetic simulations outcomes of Gagnier & Pejcha (2023), our third aim was to assess the impact of magnetic fields on the binary separation evolution timescale, on angular momentum transport within the shared envelope, on the short-term variability of mass accretion onto the binary, and on the formation of overdensities. We found that magnetic fields play little to no role in the binary separation evolution. This result may, however, be closely tied to our choice of boundary conditions for the magnetic field and may thus call for further investigation through numerical simulations on the scale of binary interaction that can accurately resolve the accretion flow and the interaction between magnetic fields and the cores. This will be explored in future works.

We found magnetic fields do not change the disk-like morphology of radial angular momentum transport found by Gagnier & Pejcha (2023). However, when accretion onto the central binary is allowed, we found that turbulent Maxwell stress has a comparable or even larger net contribution to the radial angular momentum transport than the Reynolds stress. Turbulent Maxwell stress transport of angular momentum points outward at all radii and it locally reverses the direction of angular momentum transport that is otherwise dominated by the mean-flow advection. When accretion is prevented by imposing reflecting boundary conditions, the net contribution of Maxwell stress to the radial transport of angular momentum is negligible compared to the contribution from Reynolds stress and mean-flow.

Our investigation of the intricate dynamics of nonaxisymmetric overdensities within the common envelope of binary star systems reveals their pivotal role in modulating mass accretion onto the binary components. Similarly to Gagnier & Pejcha (2023), we find mass accretion to exhibit distinct frequencies, including the 2Ω_orb frequency stemming from the quadrupolar moment of the binary’s potential and a lower-frequency ≃0.2Ω_orb modulation associated with overdensities known as “lumps” in the context of circumbinary disks. Our analysis uncovers the origins of these various overdensities: they form as a result of the intricate interplay between spiral density waves launched by the central binary and spiral density waves with a much larger pitch angle, which arise from interaction with vortices. Additionally, we unveiled the emergence of m = 1 accretion streams associated with the ≃0.1Ω_orb frequency. These streams owe their existence to the stabilizing effect of the magnetic tension from the strong toroidal field about the orbital plane, which prevents overdensities from being destroyed by turbulence and enables them to accumulate mass and eventually migrate towards the binary. Ultimately, these denser structures are destroyed by spiral waves and tidal forces in the binary’s vicinity.

Finally, Figs. 3 and 4 show a feature that looks like an “S-shaped” cavity. We saw this cavity already in our previous work (Gagnier & Pejcha 2023). This cavity results from the interplay of centrifugal force, which “pushes” material perpendicularly away from the rotation axis, and stochastic convection motion disrupting the symmetry of the resulting low-density chimney. Although our simulations do not exhibit the launching of jets or jet-like outflows, the presence of the cavity suggests the possibility of collimation of a “wobbling” jet, which could originate around one of the cores. It is exciting to speculate that “wobbling” of jets seen in many planetary nebulae might not arise due to the precession from the second core, but instead due to time-variable S-shaped cavity in the circumbinary envelope. Resolving the inner region, specifically inside the orbit, is probably necessary to determine whether outflows are actually launched. We plan to study this aspect in our future work.

This study has provided a comprehensive examination of the impact of magnetic fields on the post-dynamical inspiral phase of common envelope evolution. We have revealed the intricate interplay between magnetic energy amplification, energy reservoirs, and the dynamics of the envelope, shedding light on the mechanisms responsible for magnetic energy amplification and transfer. Our findings have also elucidated the role of magnetic fields in angular momentum transport and the formation of nonaxisymmetric overdensities within the common envelope. However, further investigation is needed to fully comprehend the role of magnetic fields in this context and their wider implications. In particular, recent circumbinary disk simulations have showed that resolving the cavity within the computational domain is essential for accurately measuring torques, orbital evolution timescales, and orbital eccentricity excitation (e.g., Tang et al. 2017; Moody et al. 2019; Muñoz et al. 2019; Muñoz & Lithwick 2020; Duffell et al. 2020; Tiede et al. 2020; Dittmann & Ryan 2021; Combi et al. 2022), and to self-consistently launch jets (e.g., Gold et al. 2014). Adopting this approach is likely even more crucial in the context of CEE because of the probable absence of a low-density cavity. Instead, we anticipate that the region inside the binary orbit exhibits highly complex gas dynamics and substantial magnetic energy amplification. We plan to explore these phenomena in future works.

Acknowledgments

We thank the anonymous referee for comments that improved this paper. We acknowledge fruitful discussions with K. Tomida, S. Takasao, G. Wafflard-Fernandez, and G. Lesur. The research of D.G. and O.P. has been supported by Horizon 2020 ERC Starting Grant ‘Cat-In-hAT’ (grant agreement no. 803158). This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254).

References

Appendix A: Horizontally integrated energy spectrum

The Fourier transform of a function f(r) = f(r, θ, φ) reads

$$\begin{array}{c}\hfill \widehat{f}(\mathit{k})=\widehat{f}(k,{\theta }_k,{\phi }_k)={\displaystyle {\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{\infty }f(\mathit{r})e^{i\mathit{k}·\mathit{r}}{r}^{2}sin\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi .}\end{array}$$(A.1)

We expand both f and the Fourier kernel e^ik ⋅ r on the spherical harmonics basis

$$\begin{array}{c}\hfill f(\mathit{r})={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{f}_m^{\ell }(r)Y_{\ell }^m(\theta ,\varphi ),}\end{array}$$(A.2)

and

$$\begin{array}{c}\hfill e^{i\mathit{k}·\mathit{r}}=4\pi {\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{i}^{\ell }{j}_{\ell }(kr)Y_{\ell }^m{(\theta ,\varphi )}^{\ast }{Y}_{\ell }^m({\theta }_k,{\varphi }_k),}\end{array}$$(A.3)

where $Y_{\mathcal{l}}^m$ is the usual scalar spherical harmonics function and j_ℓ(z) is the spherical Bessel function of the first kind of order ℓ. Injecting (A.2) and (A.3) into (A.1) yields

$$\begin{array}{c}\hfill \widehat{f}(\mathit{k})={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{F}_{\ell }^m(k)Y_{\ell }^m({\theta }_k,{\varphi }_k),}\end{array}$$(A.4)

where $F_{\mathcal{l}}^m(k)$ is the ℓ^th order spherical Bessel transform of the spherical harmonics expansion coefficients

$$\begin{array}{c}\hfill F_{\ell }^m(k)=4\pi i^{\ell }{\displaystyle {\int }_0^{\infty }{j}_{\ell }(kr)f_m^{\ell }(r)r^2\mathrm{d}r.}\end{array}$$(A.5)

The horizontally integrated (in Fourier space) energy spectrum of f(r) finally reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E_f(k)& =k^2{\displaystyle {\int }_0^{2\pi }{\int }_0^{\pi }\widehat{f}(\mathit{k})\widehat{f}{(\mathit{k})}^{\ast }sin{\theta }_k\mathrm{d}{\theta }_k\mathrm{d}{\phi }_k}\hfill \\ \hfill & =k^2{\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{F}_{\ell }^m(k)F_{\ell }^m{(k)}^{\ast }.}\hfill \end{array}\end{array}$$(A.6)

This quantity relates to the total energy through Parseval’s theorem

$$\begin{array}{c}\hfill ⟨E_f⟩={\displaystyle {\int }_{\partial V}{|f(\mathit{r})|}^2\mathrm{d}V=\frac{1}{{(2\pi )}^3}{\int }_0^{\infty }{E}_f(k)\mathrm{d}k.}\end{array}$$(A.7)

Appendix B: Magnetic energy evolution equation

The magnetic energy evolution equation (32) is derived from the dot product of B with the induction equation (4), that is

$$\begin{array}{c}\hfill \frac{1}{2}\frac{\partial (\mathit{B}·\mathit{B})}{\partial t}-\mathit{B}·(\mathit{B}·\mathbf{\nabla }\mathit{u})+\mathit{B}·\mathit{B}(\mathbf{\nabla }·\mathit{u})+\mathit{B}·(\mathit{u}·\mathbf{\nabla }\mathit{B})=0.\end{array}$$(B.1)

Using the fact that B ⋅ (u⋅∇B) = u ⋅ (∇B⋅B) and ∇B ⋅ B = B ⋅ ∇B + B × (∇×B), yields

$$\begin{array}{c}\hfill \frac{1}{2}\frac{\partial (\mathit{B}·\mathit{B})}{\partial t}-\mathit{B}·(\mathit{B}·\mathbf{\nabla }\mathit{u})+E_B\mathbf{\nabla }·\mathit{u}+\mathbf{\nabla }·\left(E_B\mathit{u}\right)=0,\end{array}$$(B.2)