Rings, Shells, and Arc Structures Around B[e] Supergiants. I. Classical Tools of Nonlinear Hydrodynamics

During their post-main-sequence evolution, massive stars (M > 8 M_☉) can undergo phases of enhanced mass loss and material ejections, which may lead to the formation of circumstellar shells and disks. One group of objects is particularly peculiar. These are the B[e] supergiants. The early findings by Zickgraf et al. (1985) of a hybrid character of the UV and optical spectra have led to the assumption of a two-component wind emanating from these objects with a classical line-driven wind in polar direction along with a much slower, cool, and dense equatorial (disk-forming) outflow in equatorial direction.

Intense infrared excess emission points to the presence of significant amounts of hot circumstellar dust (Zickgraf et al. 1986), which imprints its specific emission features on the infrared spectra of these objects (Kastner et al. 2010). The existence of a disk-shaped dusty structure has been inferred from optical linear polarization observations (Magalhaes 1992; Melgarejo et al. 2001) and has later on been reinforced for a few close-by Galactic objects by optical interferometric observations (Domiciano de Souza et al. 2007; Borges Fernandes et al. 2011; Millour et al. 2011; Cidale et al. 2012; Wheelwright et al. 2012a).

A warm and dense circumstellar disk provides an ideal environment for the formation of molecules, and rovibrational emission from hot molecular gas has been detected in the near-infrared from CO (McGregor et al. 1988; Morris et al. 1996; Miroshnichenko et al. 2005; Oksala et al. 2012; Wheelwright et al. 2012b; Kraus et al. 2013; Oksala et al. 2013; Kraus et al. 2014; Muratore et al. 2015; Sholukhova et al. 2015; Kraus et al. 2016; Kourniotis et al. 2018; Kraus et al. 2020) as well as from SiO (Kraus et al. 2015), and in the optical possibly from TiO (Zickgraf et al. 1989; Torres et al. 2012; Kraus et al. 2016; Torres et al. 2018) in a large number of objects (for a detailed description we refer the reader to the review by Kraus 2019). The detection of significant enrichment of the disk material in the isotope ¹³C, traced by intense emission of the molecule ¹³CO, reinforces that the circumstellar matter of B[e] supergiants must have been released from the stellar surface and cannot be a remnant from star formation (Kraus 2009; Liermann et al. 2010).

Several theoretical approaches have been presented to explain the formation of dense outflowing disks from B[e] supergiants. Bjorkman & Cassinelli (1993) proposed that the winds emanating from the polar regions of rapidly rotating massive stars are bent toward the equator regions where they collide forming the so-called wind-compressed disk. This analysis considered a spherical shape of the star, but the inclusion of the non-radial forces occurring on the rotationally distorted stellar surface seems to prevent the formation of a wind-compressed disk (Owocki et al. 1996).

Rapid stellar rotation facilitates another mechanism that might be considered, the rotation-induced bistability (Pelupessy et al. 2000). Due to the decrease in surface temperature from the pole to the equator caused by the rotation of the star (known as gravity darkening), the threshold temperature of ∼25,000 K for recombination of Fe iv into Fe iii might be crossed. Because Fe iii has significantly more lines suitable to drive the wind, a substantial increase in mass flux can be expected at lower temperature, that is, toward equatorial regions (Vink et al. 1999). However, the density enhancement that can be achieved by this scenario remains a factor of 10–100 below the expectations from observations. When the bistability mechanism is combined with the slow-wind solutions discovered by Curé (2004), the situation improves and a higher density contrast can be generated in the vicinity of the star (Curé et al. 2005), but for the price of an equatorial wind velocity that is 10–20 times higher than what has been inferred from observations.

The effects of gravity darkening have also been utilized in computations of the latitude-dependent ionization structure in the winds of B[e] supergiants, and Kraus (2006) has shown that, despite the fact that rapid rotation alone leads to a lower wind density in equatorial directions as was shown by Owocki et al. (1996) and Maeder & Desjacques (2001), the wind material (in particular hydrogen and elements with a similar ionization potential) can still recombine in the equatorial wind region of these luminous objects leading to a zone of neutral gas confined to the equatorial plane, in which also molecules and dust can form. The existence of such neutral (in hydrogen) zones (or disk-like structures) has been concluded from the analysis of the observed line luminosity of the [O i] lines that arise in these disks (Kraus et al. 2007). Follow-up 2D models revealed that such recombination scenarios require very high equatorial stellar mass-loss rates (Zsargó et al. 2008). In a different approach of 2D dense, viscous outflowing disks it was found that viscous heating dominates the innermost disk regions leading to extremely high temperatures within the disk midplane and to instabilities with significant waves or bumps in density and temperature (Kurfürst et al. 2018).

Many new observations have been carried out in the past years providing clearer insight with respect to the density distribution and the dynamics within the circumstellar (disk) matter. In particular, it has been found that the circumstellar material is confined in a series of rings, arcs, or spiral-arm-like structures revolving around the central object on (quasi-)Keplerian orbits, rather than being spread over a disk in the classical picture (e.g., Kraus et al. 2010; Millour et al. 2011; Aret et al. 2012; Kraus et al. 2016; Maravelias et al. 2018; Torres et al. 2018; Kraus et al. 2023). The arrangement of these rings is thereby unique for each object (Maravelias et al. 2018), and each ring can have a different density that does not necessarily follow the usual radial distribution expected in an outflow. Moreover, these rings can have gaps or inhomogeneities, and they can be either stable in time (Kraus et al. 2016, 2023) or display temporal variabilities (Maravelias et al. 2018; Torres et al. 2018) including fading (Kraus et al. 2020), complete disappearance (Liermann et al. 2014), but occasionally also a sudden appearance of a new structure (Oksala et al. 2012) possibly caused by a pileup of matter in a steadily decelerating outflow (Kraus et al. 2010). These findings indicate that the mass loss from these stars is not a smooth process, but could be related to ejection phases, possibly triggered by instabilities acting in the strongly inflated envelopes of such massive and luminous objects (Glatzel & Kiriakidis 1993; Kiriakidis et al. 1993; Glatzel & Mehren 1996).

Motivated by this great diversity of circumstellar environments of B[e] supergiants ranging from stationary density distributions in the form of rotating rings with sometimes alternating densities, or arc-like features, to decelerating equatorial outflows with sudden pileup of matter, and the deficiency of existing models to describe them, we develop in this paper new perspectives with different hydrodynamical scenarios that might help in understanding their formation ¹ and stationary structure. In particular, the formerly discussed scenarios and phenomena should be interconnected with basic, generic properties of fluid dynamics. We want to consider the following open questions from an abstract perspective:

• 1.  
  How can a stationary, ideal fluid representation (without dissipative effects) of persistent matter (stable, quasi-stationary rings, arcs, complete, or incomplete spiral structures) be constructed?
• 2.  
  How can quasi-stationary mass distributions and probably time-dependent flows (e.g., in the form of episodic mass loss of the central star) appear together?
• 3.  
  How can simplified single-fluid time-dependent and time-independent (stationary) velocity fields be constructed for such abovementioned cases, if not many detailed physical parameters are known?
• 4.  
  Do such time-dependent velocity fields exist at all for a stationary density distribution? What is their nature?

The paper is structured in the following way: in Section 2, we present the stationary 2D solution techniques but also the transition to nonstationary and 3D problems within incompressible ideal hydrodynamics (HD). In Section 3, we derive a nonlinear Schrödinger-type equation as an equivalent formulation of the compressible HD equations. In Section 4, we analyze a general 3D compressible flow on the basis of stationary stellar wind solutions. Our results are discussed and our conclusions are summarized in Section 5.

2.1. Basic Equations of Ideal HD

The scenarios described in the introduction will be treated by using the basic equations of HD. These are given by the mass continuity equation, Equation (1), and the Euler equation, Equation (2),

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}{\boldsymbol{\cdot }}\left(\rho {\boldsymbol{v}}\right)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial v}{\partial t}+\left({\boldsymbol{v}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right){\boldsymbol{v}}\right)=-{\boldsymbol{\nabla }}p+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

in which ρ is the mass density, v the gas velocity, p is the gas pressure, which may include also the radiation pressure given that it can be assumed to be isotropic ² , and g = −∇ ϕ is the gravitational acceleration of the star with the gravitational potential ϕ. Self-gravitation effects of the circumstellar matter are neglected. We focus on systems in which viscosity can be neglected, meaning that the length scales of the pressure force are much larger than the so-called deflection length, and the flows are supersonic but not highly supersonic (i.e., no shocks involved, see Frank et al. 1992).

To allow for a profound investigation of the physical aspects, we will split our analysis into two distinct physical extremes: incompressible velocity fields (Section 2.2) and irrotational potential velocity fields (Section 3). One should be aware of the fact that the incompressibility condition ∇ · v = 0 is neither valid for classical viscous disk models, nor for classical wind solutions, as these usually assume isothermy. Our incompressibility model is not isothermal, and it makes shocks less likely to occur.

2.2. Stationary 3D Incompressible Flows

An interesting case that has emerged from the observations is the (presumably) Keplerian rotating rings detected around the B[e]SG star LHA 120-S 73 (Kraus et al. 2016). These rings did not show any measurable radial motion within the observation period spanning 16 yr, justifying the assumption of a quasi-stationary model, i.e., time-dependent changes of physical quantities are small and can be neglected. Moreover, due to the preferential rotational and relaxed motion of the gas along closed orbits, the compressibility of the gas is negligible meaning that ∇ · v = 0. In this case, the mass conservation Equation (1) simplifies to

$$\begin{eqnarray}&&{\boldsymbol{v}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\rho =0,\end{eqnarray} \tag{ 3 }$$

implying that the flow is always perpendicular to the gradient of the density ( v ⊥∇ ρ), respectively, the density is constant along streamlines.

We now introduce the streaming vector (or auxiliary flow vector, Nickeler & Fahr 2005; Nickeler et al. 2006) w via

$$\begin{eqnarray}&&{\boldsymbol{w}}=\sqrt{\rho }{\boldsymbol{v}}\end{eqnarray} \tag{ 4 }$$

so that

$$\begin{eqnarray}&&{\boldsymbol{w}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\rho =0\quad \wedge \quad {\boldsymbol{\nabla }}{\boldsymbol{\cdot }}{\boldsymbol{w}}=0.\end{eqnarray} \tag{ 5 }$$

The total pressure, or as we will call this in the following, the Bernoulli pressure Π, is defined by

$$\begin{eqnarray}&&{\rm{\Pi }}=p+\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}=p+\displaystyle \frac{1}{2}{{\boldsymbol{w}}}^{2}.\end{eqnarray} \tag{ 6 }$$

Applying the identity

$$\begin{eqnarray}&&\left({\boldsymbol{w}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right){\boldsymbol{w}}={\boldsymbol{\nabla }}\left(\displaystyle \frac{1}{2}{{\boldsymbol{w}}}^{2}\right)-{\boldsymbol{w}}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{w}}\right)\end{eqnarray} \tag{ 7 }$$

to the stationary Euler equation, (Equation (2), and the definitions for the streaming vector in Equations (4)–(6), the resulting momentum equation in 3D can be written in the form

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{\rm{\Pi }}={\boldsymbol{w}}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{w}}\right)-\rho {\boldsymbol{\nabla }}\phi .\end{eqnarray} \tag{ 8 }$$

Equations (5) and (8) form the set of incompressible HD equations that need to be solved.

2.2.1. Solutions in 2D

Due to flat or disk-like structuring and concentration of the matter around the star, we assume that the material is confined in a thin sheet around the equatorial plane. In general, circumstellar disks can be assumed to be symmetric around the midplane, having their maximum density and pressure along z = 0 so that ∂/∂z = 0. Our analysis is restricted to the x–y plane and fringe effects are neglected. The investigated scenario is not intended to generate an outflowing disk in the classical sense, but to find representations for revolving rings or arcs in a quasi-stationary state. In the limit, we assume a purely azimuthal flow and use a different geometry (basically 2D in Cartesian components, i.e., (x, y) = (R, Φ)), where in contrast most works about disks use the (R, z) coordinate system (i.e., ∂/∂Φ = 0).

To satisfy the mass conservation equation, Equation (5), we define the stream function ψ = ψ(x, y), the streaming vector w = ∇ ψ × e _z , and the mass density ρ = ρ(ψ). In the general case, ∇ ψ and ∇ ϕ are not parallel to each other almost everywhere, such that the stream function ψ and the gravitational potential ϕ can be regarded as coordinates replacing x and y, and Π is consequently an explicit function of ψ and ϕ. Inserting these relations into the momentum equation, Equation (8), and expanding with respect to the basis vectors ∇ ψ and ∇ ϕ delivers

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{\rm{\Pi }}={\rm{\Delta }}\psi {\boldsymbol{\nabla }}\psi -\rho (\psi ){\boldsymbol{\nabla }}\phi ,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\Rightarrow \displaystyle \frac{\partial {\rm{\Pi }}}{\partial \psi }{\boldsymbol{\nabla }}\psi +\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \phi }{\boldsymbol{\nabla }}\phi ={\rm{\Delta }}\psi {\boldsymbol{\nabla }}\psi -\rho (\psi ){\boldsymbol{\nabla }}\phi .\end{eqnarray} \tag{ 10 }$$

From a comparison of coefficients we obtain one nonlinear Poisson-like partial differential equation for the stream function and one equation for the external gravitational potential

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \psi }={\rm{\Delta }}\psi \quad \wedge \quad \displaystyle \frac{\partial {\rm{\Pi }}}{\partial \phi }=-\rho (\psi ).\end{eqnarray} \tag{ 11 }$$

This set of relations is an analogy to results in magnetohydrostatic theory with gravity (Schindler et al. 1983), while the first quasi- or nonlinear elliptic type equation has already been found and described in HD by Stokes (1848). Formal integration of the equation on the the right-hand side of Equation (11) results in a relation for the Bernoulli pressure of

$$\begin{eqnarray}&&{\rm{\Pi }}=-\rho (\psi )\phi +{{\rm{\Pi }}}_{0}(\psi );\end{eqnarray} \tag{ 12 }$$

however, this formal integration does not solve the system completely. Inserting this Bernoulli pressure into the equation on the left-hand side of Equation (11) delivers

$$\begin{eqnarray}&&{\rm{\Delta }}\psi =-\rho ^{\prime} (\psi )\phi +{{\rm{\Pi }}}_{0}^{{\prime} }(\psi ),\end{eqnarray} \tag{ 13 }$$

where the prime denotes the derivative with respect to ψ.

To be able to solve this equation, it is necessary to know the function ρ(ψ), i.e., the density as a function of the stream function. While the stream function ψ is not known a priori, the density along each streamline label (i.e., the value of the stream function) must be constant, but can vary across streamlines. With this knowledge, we are able to calculate the 2D density structure, as we can basically choose the density function ρ = ρ(ψ) arbitrarily. The choice of the density function implies a nonlinear feedback on the solution of the nonlinear Laplace equation, Equation (13). The solution of this Laplace equation finally delivers a stream function ψ = ψ(x, y), based on which the spatial density distribution ρ = ρ(x, y) can be computed ³ .

Taking the most simple, nontrivial approach given by

$$\begin{eqnarray}&&\rho ^{\prime} (\psi )=-\lambda =\mathrm{const},\end{eqnarray} \tag{ 14 }$$

and assuming that Π₀ is constant, results in

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \psi }={\rm{\Delta }}\psi =\lambda \phi .\end{eqnarray} \tag{ 15 }$$

The relation between the density and the stream function, Equation (14), is chosen in such a way that the monotonicity of the density function, and therefore a unique relation, i.e., the bijective character of the density function, is guaranteed. We will recognize in the following that λ controls the influence of the gravitation on the geometry and dynamics of the flow: the larger λ, the larger the deviation from ordinary potential flows (see Equations (17) and (18) below).

To derive this influence we first of all have to facilitate the solution of Equation (15), by switching to Wirtinger calculus (e.g., Remmert 1991) and using the coordinates (or variables) u := x + iy and v := x − iy with i² = −1. Then the Poisson equation, Equation (15) can be written as

$$\begin{eqnarray}&&4{\partial }_{u}{\partial }_{v}\psi =-\lambda \,\displaystyle \frac{{{GM}}_{* }}{\sqrt{{uv}}},\end{eqnarray} \tag{ 16 }$$

where we inserted the definition of the gravitational potential given by $\phi =\tfrac{-{{GM}}_{* }}{R}$, in which G is the gravitational constant, M_* is the mass of the central star, and where R defined as R² = uv is the radial distance from the center of the star within the x–y-plane, with R ≥ R_*, where R_* denotes the radius of the star.

Integration of the solution of the Poisson equation, Equation (16), with respect to the u–v coordinates delivers the stream function

$$\begin{eqnarray}&&\psi =-\left(\lambda {{GM}}_{* }\right)\sqrt{{uv}}+{\psi }_{1}(u)+{\psi }_{2}(v),\end{eqnarray} \tag{ 17 }$$

where ψ₁ and ψ₂ are free functions. With this stream function, we can compute the streaming vector and the Bernoulli pressure, where the thermal pressure p can then be calculated by subtracting the kinetic pressure. The general solution, Equation (17), shows an inhomogeneous part depending on the coupling constant λ.

This simplified approach leads to a mathematical limit to a quasi-Kepler rotation. We assume that the homogeneous part of the general solution, Equation (17), i.e., the meromorphic part, is zero. Therefore the solution for ψ and the derived streaming vector can be used to calculate the (rotational) velocity of the gas around the star

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{w}} & = & -\left(\lambda {{GM}}_{* }\right){\boldsymbol{\nabla }}R{\boldsymbol{\times }}{{\boldsymbol{e}}}_{z}\\ & = & -\left(\lambda {{GM}}_{* }\right)\left[-\displaystyle \frac{x}{R}{e}_{y}+\displaystyle \frac{y}{R}{e}_{x}\right]\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\Rightarrow \,{\boldsymbol{v}}=\displaystyle \frac{{\boldsymbol{w}}}{\sqrt{\rho (\psi )}}=\displaystyle \frac{-\lambda {{GM}}_{* }}{\sqrt{{\lambda }^{2}{{GM}}_{* }R+{\rho }_{0}}}\left(\displaystyle \frac{y}{R},-\displaystyle \frac{x}{R}\right)\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\Rightarrow \,{v}^{2}=\displaystyle \frac{{{GM}}_{* }}{R+\tfrac{{\rho }_{0}}{{\lambda }^{2}{{GM}}_{* }}}=\displaystyle \frac{{{GM}}_{* }}{R}\,\displaystyle \frac{1}{1+\tfrac{{\rho }_{0}}{{\lambda }^{2}{{GM}}_{* }R}},\end{eqnarray} \tag{ 20 }$$

where we have inserted for the density function the solution of Equation (14), i.e., ρ(ψ) = −λ ψ + ρ₀ = λ² GM_* R + ρ₀. For the limit λ, ρ₀ → 0 and ρ₀/λ² → 0, or for ρ₀ = 0, the velocity is identical to the Kepler rotation. For a vanishing integration constant ρ₀, the density at the stellar surface is given by ρ(R_*) = λ² GM_* R_*. While the density of the circumstellar matter in the vicinity of the star can take small values, depending on the choice of λ, it is not trivial to generate a Keplerian rotating disk or ring detached from the stellar surface.

Other cases, e.g., for large values of ρ₀/λ², lead to strong deviation from Kepler rotation close to the star. But, even in the latter case, the velocity can approach Keplerian behavior for sufficiently large values of R. It should be noted that our approach only represents the close-by circumstellar environment, i.e., regions between the stellar surface and the first ring or onset of the disk. The linear function chosen for the density ρ(ψ) might be regarded as the lowest order part of a Taylor expansion in the region of our interest of a complex density function, and is only meant to have a prescription from a low density, e.g., close to the star (gap), to high values (ring/disk). The lack of information on the real radial density distribution from observations hampers a more precise theoretical description.

Inserting the streaming vector, Equation (18), and the density function ρ(ψ) into the equation for the pressure function, Equation (12), Equation (6) delivers

$$\begin{eqnarray}&&p=\displaystyle \frac{{{GM}}_{* }}{R}{\rho }_{0}+\displaystyle \frac{1}{2}{\lambda }^{2}{G}^{2}{M}_{* }^{2}+{{\rm{\Pi }}}_{0}.\end{eqnarray} \tag{ 21 }$$

If ρ₀ would be zero, the density would be given by λ² GM_* R, the pressure gradient (and therefore the pressure force) would vanish, and the motion of the gas would be purely ballistic (Keplerian).

2.2.2. On the Existence of Π(ψ, ϕ) for 1D Equilibria

In Section 2.2.1, we showed for two-dimensional equilibria that the pressure Π can be expressed as a function of the scalar fields ψ and ϕ with

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \psi }={\rm{\Delta }}\psi := -{{\rm{\Omega }}}_{w}(\psi ,\phi )\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \phi }={\rm{\Delta }}\phi := -\rho (\psi ),\end{eqnarray} \tag{ 23 }$$

where Ω_w is the vortex strength (or vorticity in 2D) associated with the flow vector w, which is defined by Ω_w $:=$ ∇ × w and reduces in our case to Ω_w = Ω_w e _z .

In the one-dimensional case, e.g., ∂/∂Φ = 0, where (x, y) ≡ (R, Φ) (with R and Φ being the radial coordinate and the azimuthal angle, respectively), it is not obvious at first that there is an incompressible stationary solution, characterized by Π(R), ψ(R) and ϕ(R), the vorticity Ω_w (R), and the mass density ρ(R).

A possible proof has been proposed by Hornig (1996, private communication) and presented by Fleischer (1996), which we recapitulate for clarification: We consider the two-dimensional space ${{\mathbb{R}}}^{2}$ with the coordinates (ψ, ϕ), in which the solution ψ(R), ϕ(R) is a curve L, parameterized by R. On L the vector field V = V_ψ (ψ, ϕ) e _ψ + V_ϕ (ψ, ϕ) e _ϕ is given by

$$\begin{eqnarray}&&{V}_{\psi }(\psi (R),\phi (R))=-{{\rm{\Omega }}}_{w}\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&{V}_{\phi }(\psi (R),\phi (R))=-\rho .\end{eqnarray} \tag{ 25 }$$

We are looking for a potential Π on ${{\mathbb{R}}}^{2}$ such that V = ∇Π∣_L .

We assume that there is a ε₀-hose around L that does not overlap anywhere. Then there are local coordinates (s, ε), where s is the arc length along L and ε is the distance from L. These are orthogonal to L, e _s · e _ε ∣_L = 0. There is then a decomposition of V (s) = V_s (s) e _s + V_ε (s) e _ε on L, and the potential Π, with Π ≡ 0 outside the tube, and

$$\begin{eqnarray}&&{\rm{\Pi }}(\varepsilon ,s)=\left[{\int }_{0}^{s}\,{V}_{s}(\tilde{s})d\tilde{s}+\varepsilon {V}_{\varepsilon }(s)\right]\exp \left(-\,\displaystyle \frac{{\varepsilon }^{2}}{{\varepsilon }_{0}^{2}-{\varepsilon }^{2}}\right)\end{eqnarray} \tag{ 26 }$$

for ε ≤ ε₀ fulfills the requirement V = ∇Π∣_L inside the tube.

Thus, there is a pressure function Π(ψ, ϕ) whose partial derivatives according to Equations (22) and (23) agree with the vorticity Ω_w and the mass density ρ as functions of ψ and ϕ.

2.2.3. Preliminary Discussion: Do Time-dependent Incompressible Flows Exist in 2D or 3D?

Even if a great variety of geometrical structures with regard to flow patterns can be generated in 2D with the method described in Section 2.2.1, we know that even quasi-stationary structures have to undergo a formation phase that can be quite eruptive and therefore strongly time dependent. This leads to the question of whether quasi-stationary structures can be preserved in case of a nonlinear time-dependent change of the system or relax, respectively, blow-up. Or, equivalently, is it at least possible to extend the stationary models via slight amplitude modulation, i.e., the assumption that stationary fields can be multiplied with a specific purely time-dependent function, to deliver a time-dependent solution including the stationary solution? To answer this question we focus on the dynamical behavior in the nonlinear case. We additionally assume that the system reacts instantaneously on nonlocal perturbations.

The slight amplitude modulation is equivalent to posing the question if a time separability from the spatial components guarantees regular solutions within a small, finite-time interval. The regularity would at least make it plausible that the solutions are nonlinearly stable. The separability of solutions of evolutionary equations with respect to time is possible for many examples in continuum and quantum physics (see, e.g., Galaktionov & Svirshchevskii 2007). The general solutions of such kind of equations can then be constructed by linear superposition of multiple summands constructed by such a separation ansatz, or by nonlinear superposition of a similar set of summands.

We assume a time separability and take the general incompressible 3D representation into account. Considering the Euler equation for vanishing time derivative ∂ρ/∂t ≈ 0 and external force, and using the definition of the stream vector (Equations (4) and (5)), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{\nabla }}p & = & -\rho \left({\boldsymbol{v}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right){\boldsymbol{v}}-\rho \displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}\\ \iff \quad {\boldsymbol{\nabla }}{\rm{\Pi }} & = & {\boldsymbol{w}}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{w}}\right)-\sqrt{\rho }\,\displaystyle \frac{\partial {\boldsymbol{w}}}{\partial t}.\end{array}\end{eqnarray} \tag{ 27 }$$

With w = w₁(t) w₀(x), and analogously for the total pressure Π, we can rewrite the second form of the Euler Equation (27) as

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{1}(t){\boldsymbol{\nabla }}{{\rm{\Pi }}}_{0}({\boldsymbol{x}})={w}_{1}^{2}\,{{\boldsymbol{w}}}_{0}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{w}}}_{0}\right)-{\dot{w}}_{1}\,\sqrt{\rho }\,{{\boldsymbol{w}}}_{0},\end{eqnarray} \tag{ 28 }$$

where the dot in ${\dot{w}}_{1}$ denotes the derivative with respect to time. Assuming that ∃( w ₀, Π₀, ρ) with

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{{\rm{\Pi }}}_{0}={{\boldsymbol{w}}}_{0}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{w}}}_{0}\right)-\sqrt{\rho }\,{{\boldsymbol{w}}}_{0}\end{eqnarray} \tag{ 29 }$$

delivers

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{1}={w}_{1}^{2}\propto {\dot{w}}_{1}\quad \Rightarrow \quad {w}_{1}=\displaystyle \frac{{w}_{10}}{{t}_{0}-t},\end{eqnarray} \tag{ 30 }$$

with the integration constant w₁₀. For w₁₀ > 0, t₀ > 0, and t₀ > t this solution constitutes a blow-up solution or finite-time singularity for t → t₀, whereas in the limit t → − ∞ the velocity and pressure decay. We would like to emphasize that Equation (29) does not represent an equilibrium solution, but poses an additional constraint for the stationary part of the flow, respectively the steady-state flow pattern.

If we take also gravity into account for the above calculations then the time-independent form of the momentum equation is

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{\nabla }}{{\rm{\Pi }}}_{0} & & ={{\boldsymbol{w}}}_{0}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{w}}}_{0}\right)-\sqrt{\rho }\,{{\boldsymbol{w}}}_{0}-\displaystyle \frac{{{GM}}_{* }\rho }{{r}^{2}}\,{{\boldsymbol{e}}}_{r}\\ \Rightarrow \quad {{\rm{\Pi }}}_{1} & & ={w}_{1}^{2}={\dot{w}}_{1}\propto {M}_{* }(t),\end{array}\end{eqnarray} \tag{ 31 }$$

where e _r is the unit vector directing radially outward from the center of gravity, and M_*(t) is the stellar mass changing over time. For t₀ < 0 and t = 0 the relation ${\dot{w}}_{1}={w}_{10}/{t}_{0}^{2}=1$ fulfills the condition that M_*(t = 0) = M_*. Then it follows that for t > 0 and t → ∞ the mass, pressure, and velocity decrease, whereas they increase for t₀ > 0, t₀ > t, and t → t₀. While the former can be interpreted with mass loss, the latter would imply mass accretion. It should be noted that both scenarios only hold for a meaningful (small) time interval, as the temporal variation of the velocity field is given by a (pure) dilatation.

The assumptions leading to Equations (27)–(31) reflect a nonlinear perturbation of a system, whose original stationary (geometrical) structure should be preserved, whereas the amplitudes of the stationary field components should vary only slowly. The calculations show that the time-independent field components can in fact no longer satisfy any stationary equation, and that even the time-dependent amplitudes diverge, i.e., develop finite-time singularities, or decay completely. Although the dynamics are limited to incompressible flows in these considerations, the behavior of the time-dependent amplitudes implies an eruptive mass loss (or mass gain) of the star. This means that within a short period of time the stellar mass loss (or gain) would be considerably enhanced when the system is forced by nonlinear variations (perturbations), which can change the character of the system, driving it away from a relaxed incompressible flow. Our analysis shows that time-dependent incompressible flows around a stationary state cannot exist but they lead either to a completely decaying flow or to a blow-up.

In contrast to the studies in Section 2.2, the time-dependent change in the velocity field is now included. Also, we drop the restriction to incompressible flows. Instead, we assume that the flow is irrotational. Moreover, we replace the stream function with a scalar potential, which allows us to easily implement the time dependence of the velocity field.

In contrast to the stream function model used in Section 2.2, reflecting the isocontours of the density distribution in a relaxed flow, we now utilize a scalar velocity potential. The surfaces of constant potential represent a family of surfaces, which has an affinity to the family of surfaces of constant radial coordinate. This would favor a radial outflow, instead of the more azimuthal flows described by the stream function model in the previous section. In case of a time dependency, the potential describes expanding (or eventually shrinking) surfaces, to which the velocity vector is perpendicular, thus pointing in a certain sense radially outward (or inward). The density, decreasing outwards (eventually inward), should be therefore, at least locally, diffeomorphic with respect to the potential, as will be introduced and explained in the following.

First of all, maintaining the condition of a stationary density distribution, the equation of mass conservation reduces again to

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{\boldsymbol{\cdot }}\left(\rho {\boldsymbol{v}}\right)=0.\end{eqnarray} \tag{ 32 }$$

The Euler equation, Equation (2), can be rewritten, using the identity equation, Equation (7), as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{\boldsymbol{\nabla }}\displaystyle \frac{{{\boldsymbol{v}}}^{2}}{2}-{\boldsymbol{v}}{\boldsymbol{\times }}\left({\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{v}}\right)=-\displaystyle \frac{{\boldsymbol{\nabla }}p}{\rho }-{\boldsymbol{\nabla }}\phi .\end{eqnarray} \tag{ 33 }$$

Introducing a 3D velocity potential v := ∇ φ, leading to ∇ × v = 0, and utilizing a barotropic law p = p(ρ), Equations (32) and (33) take the following form:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\varphi {\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\rho =-\rho {\rm{\Delta }}\varphi ,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{\nabla }}\left(\int \frac{{dp}}{\rho (p)}+\frac{{\left({\boldsymbol{\nabla }}\varphi \right)}^{2}}{2}+\frac{\partial \varphi }{\partial t}+\phi \right)={\bf{0}}\,\\ \iff \qquad \,\,\int \frac{{dp}}{\rho (p)}+\frac{{\left({\boldsymbol{\nabla }}\varphi \right)}^{2}}{2}+\frac{\partial \varphi }{\partial t}+\phi =0.\end{array}\end{eqnarray} \tag{ 35 }$$

The latter, Equation (35), is the classical Bernoulli equation for compressible, unsteady flows without viscosity. Without loss of generality, we have incorporated a possible time-dependent integration constant into the function φ. This equation needs further integration, as the mass conservation equation, Equation (34), must also be fulfilled. The study of possible temporal evolutions of the Euler equation, Equation (33), together with the mass continuity equation, Equation (34), is the subject of Section(3.1).

In analogy to the incompressibility representation (ρ = ρ(ψ)), we introduce the affinity between isosurfaces of the velocity potential and density, namely, ρ = ρ(φ), where $\rho ^{\prime} =\rho ^{\prime} (\varphi )=d\rho /d\varphi$. Adding an explicit time dependence to the mass continuity equation, Equation (34), results in the time-dependent continuity equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial \varphi }{\partial t}+{\left({\boldsymbol{\nabla }}\varphi \right)}^{2}+\displaystyle \frac{\rho }{\rho ^{\prime} }\,{\rm{\Delta }}\varphi =0.\end{eqnarray} \tag{ 36 }$$

Using this equation, we substitute either the quadratic derivative of the potential φ or its time derivative in the Euler equation, Equation (35), and obtain

$$\begin{eqnarray}&&G(\varphi )+\displaystyle \frac{1}{2}\,\displaystyle \frac{\partial \varphi }{\partial t}-\displaystyle \frac{\rho }{2\rho ^{\prime} }\,{\rm{\Delta }}\varphi +\phi =0,\end{eqnarray} \tag{ 37 }$$

respectively,

$$\begin{eqnarray}&&G(\varphi )-\displaystyle \frac{{\left({\boldsymbol{\nabla }}\varphi \right)}^{2}}{2}-\displaystyle \frac{\rho }{\rho ^{\prime} }\,{\rm{\Delta }}\varphi +\phi =0\end{eqnarray} \tag{ 38 }$$

with

$$\begin{eqnarray}&&G(\rho (\varphi ))=\int \displaystyle \frac{{dp}}{\rho (p)}=\int \displaystyle \frac{{dp}(\rho )}{d\rho }\,\displaystyle \frac{d\rho }{d\varphi }\,\displaystyle \frac{d\varphi }{\rho (\varphi )}.\end{eqnarray} \tag{ 39 }$$

The function G can be identified as the specific enthalpy.

Equation (37) is a nonlinear diffusion equation, similar to the nonlinear Schrödinger equation, and Equation (38) is a convection–diffusion-type equation.

In the following, we first discuss a special case with a specific time dependency where only φ is time dependent and ρ is stationary and depends only on the spatial parts of φ (Section 3.1), whereas in Section 3.2 we return to the general case of a density depending explicitly on the time-dependent velocity potential.

3.1. Separable Time Dependence and Neglection of Gravity

We want to search for the general time dependence of (circumstellar) flows in order to find out whether instabilities or collapse processes necessarily occur, even if the system stipulates a stationary density at far distances from the central star as observations suggest. In particular, we are interested in an expansion of a solution around an equilibrium state, and this can be naturally achieved by the separation of the temporal and spatial parts of the corresponding fields. Therefore, we assume a separability of time dependency for the pressure and velocity potential, consider a stationary density, depending on the spatial part of the velocity potential, and neglect gravity, ⁴ i.e., p = p₀(t)p₁(φ₁), φ = φ₀(t)φ₁(x, y, z), and ρ = ρ(φ₁). The chosen time separability for the pressure enables a barotropic law at any time. Inserting this separability ansatz into the conservation of mass equation for stationary density, Equation (34), we obtain the relation

$$\begin{eqnarray}&&{\left({\boldsymbol{\nabla }}{\varphi }_{1}\right)}^{2}=-\displaystyle \frac{\rho ({\varphi }_{1})}{\rho ^{\prime} ({\varphi }_{1})}{\rm{\Delta }}{\varphi }_{1}.\end{eqnarray} \tag{ 40 }$$

With the 3D velocity potential v = ∇ φ, and neglecting gravity, Equation (33) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{\boldsymbol{\nabla }}\displaystyle \frac{{{\boldsymbol{v}}}^{2}}{2}=-\displaystyle \frac{{\boldsymbol{\nabla }}p}{\rho }.\end{eqnarray} \tag{ 41 }$$

Inserting the separation ansatz for the pressure and the velocity potential leads to

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}({\dot{\varphi }}_{0}(t){\varphi }_{1})+{\boldsymbol{\nabla }}\left[\displaystyle \frac{{\varphi }_{0}^{2}(t)}{2}{\left({\boldsymbol{\nabla }}{\varphi }_{1}\right)}^{2}\right]=-\displaystyle \frac{{\boldsymbol{\nabla }}\left({p}_{0}(t){p}_{1}({\varphi }_{1})\right)}{\rho ({\varphi }_{1})}.\end{eqnarray} \tag{ 42 }$$

The right-hand side can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{\nabla }}\left({p}_{0}(t){p}_{1}({\varphi }_{1})\right)}{\rho ({\varphi }_{1})}={\boldsymbol{\nabla }}\left[{p}_{0}(t)\int \displaystyle \frac{{p}_{1}^{\prime} ({\varphi }_{1})d{\varphi }_{1}}{\rho ({\varphi }_{1})}\right].\end{eqnarray} \tag{ 43 }$$

With this, Equation (42) can be formally integrated, and if we replace the term ${({\boldsymbol{\nabla }}{\varphi }_{1})}^{2}$ by Equation (40), we obtain

$$\begin{eqnarray}&&{p}_{0}(t)\int \displaystyle \frac{{p}_{1}^{\prime} ({\varphi }_{1})d{\varphi }_{1}}{\rho ({\varphi }_{1})}-\displaystyle \frac{{\varphi }_{0}^{2}(t)}{2}\,\displaystyle \frac{\rho ({\varphi }_{1})}{\rho ^{\prime} ({\varphi }_{1})}\,{\rm{\Delta }}{\varphi }_{1}+{\dot{\varphi }}_{0}(t){\varphi }_{1}=0.\end{eqnarray} \tag{ 44 }$$

We define

$$\begin{eqnarray}&&\tilde{F}:= \int \displaystyle \frac{{p}_{1}^{\prime} ({\varphi }_{1})d{\varphi }_{1}}{\rho ({\varphi }_{1})}.\end{eqnarray} \tag{ 45 }$$

As one can see by inspection of Equation (39), $\tilde{F}$ is up to a factor p₀(t) identical to the specific enthalpy G. Inserting Equation (45) into Equation (44), we get

$$\begin{eqnarray}&&{p}_{0}-{\varphi }_{0}^{2}\,\displaystyle \frac{\rho }{\rho ^{\prime} }\,\displaystyle \frac{{\rm{\Delta }}{\varphi }_{1}}{2\tilde{F}}+\displaystyle \frac{{\dot{\varphi }}_{0}{\varphi }_{1}}{\tilde{F}}=0\,\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&\iff \qquad \displaystyle \frac{{p}_{0}}{{\dot{\varphi }}_{0}}-\displaystyle \frac{{\varphi }_{0}^{2}}{{\dot{\varphi }}_{0}}\,\displaystyle \frac{\rho }{\rho ^{\prime} }\,\displaystyle \frac{{\rm{\Delta }}{\varphi }_{1}}{2\tilde{F}}+\displaystyle \frac{{\varphi }_{1}}{\tilde{F}}=0.\end{eqnarray} \tag{ 47 }$$

Equation (47) is only valid for ${\dot{\varphi }}_{0}\ne 0$. To achieve the separability of this equation into spatial and time-dependent parts, we factorize the first terms in Equation (47) and recognize that it is indispensable that the term $\tfrac{{p}_{0}}{{\varphi }_{0}^{2}}$ must be constant. We label this constant as c₁ and obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\varphi }_{0}^{2}}{\dot{{\varphi }_{0}}}\left({c}_{1}-\displaystyle \frac{\rho }{\rho ^{\prime} }\,\displaystyle \frac{{\rm{\Delta }}{\varphi }_{1}}{2\tilde{F}}\right)+\displaystyle \frac{{\varphi }_{1}}{\tilde{F}}=0.\end{eqnarray} \tag{ 48 }$$

Equation (48) can only be solved if $\tfrac{{\varphi }_{0}^{2}}{{\dot{\varphi }}_{0}}={\varphi }_{00}={\rm{const}}$. The integral of this differential equation is

$$\begin{eqnarray}&&{\varphi }_{0}=\displaystyle \frac{{\varphi }_{00}}{{t}_{0}-t},\end{eqnarray} \tag{ 49 }$$

delivering a finite-time singularity for the time-dependent part of the velocity potential. The constant φ₀₀ describes the amplitude of the velocity potential imprinted onto the spatial part of the velocity potential for the time t = 0.

For the spatial part, we find

$$\begin{eqnarray}&&2\displaystyle \frac{{\varphi }_{1}}{{\varphi }_{00}}\displaystyle \frac{\rho ^{\prime} }{\rho }+2\,{c}_{1}\tilde{F}\displaystyle \frac{\rho ^{\prime} }{\rho }={\rm{\Delta }}{\varphi }_{1}.\end{eqnarray} \tag{ 50 }$$

Equation (50) is basically a nonlinear elliptic partial differential equation whose solution reflects the spatial part of the scalar velocity potential φ₁ and thus determines the structure of the entire time-dependent flow. To solve this equation requires the knowledge of the specific enthalpy $\tilde{F}$, respectively, G, which results from the specification of the density function ρ(φ₁) along with pressure dependency p₁(ρ) or explicitly p₁(φ₁).

A detailed analysis of this equation is beyond the scope of the current paper, however, we wish to draw attention to the fact that, if the time separability is requested along with a time independent density distribution (stationarity), then according to Equation (49) for φ₀₀ > 0, t₀ > 0 and t₀ > t the system inevitably develops into a finite-time singularity for t → t₀, whereas in the limit t → −∞ the velocity and pressure decay.

For a simplified case we can find a subspace. For example, for the case ${\dot{\varphi }}_{0}=0$, Equation (46) takes the following form:

$$\begin{eqnarray}&&{p}_{0}-{\varphi }_{0}^{2}\displaystyle \frac{\rho }{\rho ^{\prime} }\,\displaystyle \frac{{\rm{\Delta }}{\varphi }_{1}}{2\tilde{F}}=0.\end{eqnarray} \tag{ 51 }$$

This equation can only be fulfilled for ${p}_{0}={\rm{const}}$. In this case, the equation takes the form of a quasi-linear partial differential equation

$$\begin{eqnarray}&&{\rm{\Delta }}{\varphi }_{1}=2\,\displaystyle \frac{\rho ^{\prime} ({\varphi }_{1})}{\rho ({\varphi }_{1})}\displaystyle \frac{{p}_{0}}{{\varphi }_{0}^{2}}\,\tilde{F}({\varphi }_{1}),\end{eqnarray} \tag{ 52 }$$

which can also be derived directly from the nonlinear Schrödinger equation, Equation (37). Equation (52) is a nonlinear Poisson equation for which exact analytical solutions are known in 2D.

3.2. General Time-dependent Approach and Generalized Cole–Hopf Transformation

We return now to the system of nonlinear diffusion equations, Equations (36)–(39), which we want to reformulate in a more compact form. For this, we use a generalized form of the Cole–Hopf transformation. The basic form of this transformation type has been invented by Hopf (1950) and Cole (1951).

In the case of nonlinear diffusion equations, the spatial part of the differential operators can have the form

$$\begin{eqnarray}&&a{\rm{\Delta }}\varphi +b{\left({\boldsymbol{\nabla }}\varphi \right)}^{2}=c\quad \wedge \quad a,b,c\,{\rm{functions}}.\end{eqnarray} \tag{ 53 }$$

Let us assume φ = F(Λ), we can rewrite Equation (53) as

$$\begin{eqnarray}&&a\left[F^{\prime\prime} ({\rm{\Lambda }}){\left({\boldsymbol{\nabla }}{\rm{\Lambda }}\right)}^{2}+F^{\prime} ({\rm{\Lambda }}){\rm{\Delta }}{\rm{\Lambda }}\right]+b{\left(F^{\prime} ({\rm{\Lambda }})\right)}^{2}{\left({\boldsymbol{\nabla }}{\rm{\Lambda }}\right)}^{2}=c,\end{eqnarray} \tag{ 54 }$$

where the primes at F denote derivatives with respect to Λ. The transformation φ = F(Λ) should be now done in such a way that terms with (∇Λ)² are eliminated. This demand can be formulated as

$$\begin{eqnarray}&&{aF}^{\prime\prime} ({\rm{\Lambda }})+b{\left(F^{\prime} ({\rm{\Lambda }})\right)}^{2}=0.\end{eqnarray} \tag{ 55 }$$

To solve Equation (55), it is necessary that the function a/b only explicitly depends on Λ. The transformation F can be derived by integration of the condition Equation (55), which has the general solution

$$\begin{eqnarray}&&F=\int \displaystyle \frac{d{\rm{\Lambda }}}{\int \tfrac{b}{a}\,d{\rm{\Lambda }}}\,\iff \,{\rm{\Lambda }}=\int \exp \left(\int \displaystyle \frac{b}{a}\,d\varphi \right)d\varphi .\end{eqnarray} \tag{ 56 }$$

With the condition, Equation (55), Equation (54) reduces to

$$\begin{eqnarray}&&{aF}^{\prime} ({\rm{\Lambda }}){\rm{\Delta }}{\rm{\Lambda }}=c.\end{eqnarray} \tag{ 57 }$$

From a comparison of coefficients between Equations (36) and (53), we can derive the functions a, b, and c

$$\begin{eqnarray}&&a=\displaystyle \frac{\rho }{\rho ^{\prime} }\qquad {\rm{and}}\qquad b=1\qquad {\rm{and}}\qquad c=-\displaystyle \frac{\partial \varphi }{\partial t}.\end{eqnarray} \tag{ 58 }$$

Inserting the general solution, Equation (56), along with the functions a, b, and c into Equation (57) leads to

$$\begin{eqnarray}\begin{array}{rcl}\varphi ({\rm{\Lambda }}) & = & \displaystyle \int {\left[\displaystyle \int \displaystyle \frac{\rho \,d{\rm{\Lambda }}}{\rho ^{\prime} }\right]}^{-1}d{\rm{\Lambda }}\iff {\rm{\Lambda }}(\varphi )=\displaystyle \int \rho \,d\varphi \\ \Rightarrow {\rm{\Delta }}{\rm{\Lambda }} & = & -\displaystyle \frac{\partial {\rm{\Lambda }}}{\partial t}\,\displaystyle \frac{\rho ^{\prime} }{\rho }.\end{array}\end{eqnarray} \tag{ 59 }$$

The last equation of Equation (59) is a nonlinear diffusion equation and is the Cole–Hopf transformed mass continuity equation. For this kind of nonlinear diffusion or Schrödinger equations, Equations (59) and (37), it is known that they can have blow-up solutions (see, e.g., Galaktionov & Vázquez 2002).

Next, we transform Equation (38). Due to the mathematical similarity of Equations (38) and (53) we can use the transformation equation, Equation (56), for

$$\begin{eqnarray}&&a=-\displaystyle \frac{\rho }{\rho ^{\prime} }\qquad {\rm{and}}\qquad b=-\displaystyle \frac{1}{2}.\end{eqnarray} \tag{ 60 }$$

This leads to

$$\begin{eqnarray}&&\begin{array}{l}G\left(F\left({\rm{\Lambda }}\right)\right)-\displaystyle \frac{\rho }{\rho ^{\prime} }\,F^{\prime} ({\rm{\Lambda }}){\rm{\Delta }}{\rm{\Lambda }}+\phi =0\\ \Rightarrow \quad G\left(F\left({\rm{\Lambda }}\right)\right)-\displaystyle \frac{\rho }{\rho ^{\prime} }{\left[\displaystyle \int \,\displaystyle \frac{\rho ^{\prime} }{2\rho }\,d{\rm{\Lambda }}\right]}^{-1}\,{\rm{\Delta }}{\rm{\Lambda }}+\phi =0\\ \Rightarrow \quad G\left(F\left({\rm{\Lambda }}\right)\right)-\displaystyle \frac{\rho }{\sqrt{\rho }\rho ^{\prime} }{\rm{\Delta }}{\rm{\Lambda }}+\phi =0\,\\ \Rightarrow \quad G\left(F\left({\rm{\Lambda }}\right)\right)-{\left[\displaystyle \frac{d\rho }{d{\rm{\Lambda }}}\right]}^{-1}{\rm{\Delta }}{\rm{\Lambda }}+\phi =0.\end{array}\end{eqnarray} \tag{ 61 }$$

The advantage of carrying out the transformation is that Equation (61) does not contain time derivatives and quadratic terms anymore, and that it is again a nonlinear Poisson equation, which can be solved with standard procedures. The only problematic and nontrivial issue is that ρ(φ) must be chosen adequately in order to find solutions to the equation, based on which it will be possible to obtain the velocity field, pressure, and density distribution.

Now we drop the condition for incompressible or irrotational flows and treat the hydrodynamical problem of nonlinear and nonlocal instabilities from a more general point of view. It is possible to implement the gravitation of the star into the pressure in case that either p can be considered as a function of ρ (barotropic fluid), or ρ as a function of ϕ, although in the former case, it will then be necessary to slightly modify the equation of motion. Without significant loss of generality, we therefore use only one pressure gradient force and write the Euler equation for the stellar wind in the following form:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}p=-\rho \left({\boldsymbol{v}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right){\boldsymbol{v}}-\rho \displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}.\end{eqnarray} \tag{ 62 }$$

For this, the equilibrium solution, which one may also consider as the ground state of the system, is

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{p}_{0}=-{\rho }_{0}\left({{\boldsymbol{v}}}_{0}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right){{\boldsymbol{v}}}_{0}\quad \wedge \quad {\boldsymbol{\nabla }}{\boldsymbol{\cdot }}\left({\rho }_{0}{{\boldsymbol{v}}}_{0}\right)=0.\end{eqnarray} \tag{ 63 }$$

If we assume the separability for the pressure p = p₁(t)p₀( x ) and analogously also for the density ρ = ρ₁(t)ρ₀( x ) and for the velocity field v = v₁(t) v ₀( x ), then we will show in the following that for some choices of constraints there can exist regular solutions, and other choices can lead to a blow-up of the stationary nonlinear spatial solution.

The separation implies that locally, i.e., inside some small subset of a time interval, the fluid variables are not to be considered as fast evolving, around a known stationary wind solution. The mass continuity equation can then be written as

$$\begin{eqnarray}&&{\rho }_{0}{\dot{\rho }}_{1}+{\rho }_{1}{v}_{1}\,{\boldsymbol{\nabla }}{\boldsymbol{\cdot }}\left({\rho }_{0}{{\boldsymbol{v}}}_{0}\right)=0,\end{eqnarray} \tag{ 64 }$$

and with the conditions of Equation (63) it follows that

$$\begin{eqnarray}&&{\dot{\rho }}_{1}=0\qquad \Rightarrow \qquad {\rho }_{1}=1\end{eqnarray} \tag{ 65 }$$

without loss of generality.

The Euler equation, Equation (62), together with the initial condition Equation (63) can be expressed in the form given in Equation (66). Taking first the divergence operator on both sides (∇ · ), and second the curl operator (∇ × ), one receives

$$\begin{eqnarray}&&{p}_{1}\,{\boldsymbol{\nabla }}{p}_{0}={v}_{1}^{2}\,{\boldsymbol{\nabla }}{p}_{0}-{\rho }_{0}{\dot{v}}_{1}{{\boldsymbol{v}}}_{0}\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&\Rightarrow \quad {\rm{\Delta }}{p}_{0}=0\quad \wedge \quad {\boldsymbol{\nabla }}{\boldsymbol{\times }}\left({\rho }_{0}{{\boldsymbol{v}}}_{0}\right)={\bf{0}}\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&\begin{array}{l}\Rightarrow \quad \exists \,{\varphi }_{0},{\rm{\Delta }}{\varphi }_{0}=0,{v}_{00}={\rm{const.}},\,{\rm{with}}\\ {p}_{0}={v}_{00}{\varphi }_{0}\,\wedge {v}_{00}{\boldsymbol{\nabla }}{\varphi }_{0}\equiv {v}_{00}{\rho }_{0}{{\boldsymbol{v}}}_{0}={\boldsymbol{\nabla }}{p}_{0}\end{array}\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&\Rightarrow \quad {p}_{1}-{v}_{1}^{2}=-\displaystyle \frac{1}{{v}_{00}}{\dot{v}}_{1}.\end{eqnarray} \tag{ 69 }$$

The condition v₀₀ < 0 is subject to the physical approach that the direction of the velocity should be antiparallel to the pressure gradient in the relation between the mass-current vector ρ₀ v ₀, and the pressure gradient Equation (68). The assumption that v₀₀ < 0 reflects the decrease of the pressure from upwind to downwind direction, which usually drives a stellar wind or emulates decretion.

Equation (69) represents an ordinary differential equation for the time-dependent dynamics of the flow, based on the equilibrium solution Equation (63). Already under simplified assumptions, like ${p}_{1}={\rm{const.}}$, we can find solutions with finite-time singularities (blow-up solutions), but these solutions originate from a finite-time singularity and converge toward the equilibrium solution. In contrast to the decretion solution, the accretion solution would originate from the equilibrium solution and converge to a finite-time singularity. In the case of p₁ = 0, it can clearly be seen and easily calculated that v₁ ∝ 1/(t₀ − t). This solution must be excluded, as this would imply that the time perturbation reduces the complete pressure to zero for all times, implying a pressureless fluid, i.e., a completely force-free flow.

For ${p}_{1}={\rm{const.}}\ne 0$, we calculate the formal solution:

$$\begin{eqnarray}&&\begin{array}{l}t-{t}_{0}={\int }_{{v}_{10}}^{{v}_{1}}\displaystyle \frac{-\tfrac{1}{{v}_{00}}\,d{\tilde{v}}_{1}}{{p}_{1}-{\tilde{v}}_{1}^{2}}\\ \,=\displaystyle \frac{1}{2| {v}_{00}| \sqrt{{p}_{1}}}\mathrm{ln}\displaystyle \frac{{\tilde{v}}_{1}+\sqrt{{p}_{1}}}{{\tilde{v}}_{1}-\sqrt{{p}_{1}}}{| }_{{v}_{10}}^{{v}_{1}}\end{array}\end{eqnarray} \tag{ 70 }$$

leading to

$$\begin{eqnarray}&&| {v}_{00}| \sqrt{{p}_{1}}(t-{t}_{0})=\displaystyle \frac{1}{2}\left(\mathrm{ln}\displaystyle \frac{{v}_{1}+\sqrt{{p}_{1}}}{{v}_{1}-\sqrt{{p}_{1}}}-\mathrm{ln}\displaystyle \frac{{v}_{10}+\sqrt{{p}_{1}}}{{v}_{10}-\sqrt{{p}_{1}}}\right).\end{eqnarray} \tag{ 71 }$$

For t > t₀ and ${v}_{1},{v}_{10}\gt \sqrt{{p}_{1}}$ Equation (71) can be reformulated as

$$\begin{eqnarray}&&{\rm{arcoth}}\left(\displaystyle \frac{{v}_{1}}{\sqrt{{p}_{1}}}\right)=| {v}_{00}| \sqrt{{p}_{1}}(t-{t}_{0})+{C}^{2}\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&\Rightarrow \quad {v}_{1}(t)=\sqrt{{p}_{1}}\coth \left(| {v}_{00}| \sqrt{{p}_{1}}(t-{t}_{0})+{C}^{2}\right)\end{eqnarray} \tag{ 73 }$$

with

$$\begin{eqnarray}&&{C}^{2}={\rm{ar}}\coth \left(\displaystyle \frac{{v}_{10}}{\sqrt{{p}_{1}}}\right).\end{eqnarray} \tag{ 74 }$$

Then, the solution of Equation (73) implies that the velocity decays with time and relaxes into an equilibrium state.

If we assume that p₁ is not constant anymore but p₁ = p₁(t) can be imposed, then Equation (69) turns into an ordinary differential equation of Riccati type.

As the problem is nonlinear anyway, the function p₁ can be regarded, especially in this situation of a dependency of only one variable, namely, t, as a function of v₁, i.e., p₁ = p₁(v₁). For certain constraints, we are able to show that equilibria could converge to blow-up solutions. The assumption that p₁(v₁) ≥ 0 guarantees that the pressure is always positive. For a continuous outflow, we must also guarantee that ${\dot{v}}_{1}$ is always larger than zero if we exclude oscillatory phases of acceleration and deceleration. These conditions are fulfilled by ${p}_{1}\mbox{--}{v}_{1}^{2}\geqslant 0$ (see Equation (69)). Thus, it is further necessary that p₁ is not only positive, and a sufficient criterion to guarantee the validity of ${p}_{1}\mbox{--}{v}_{1}^{2}\geqslant 0$ is that there exists a p₁₀ > 1 and ${p}_{1}={p}_{10}{v}_{1}^{2}$ such that

$$\begin{eqnarray}&&{p}_{10}{v}_{1}^{2}-{v}_{1}^{2}=\displaystyle \frac{1}{| {v}_{00}| }{\dot{v}}_{1}.\end{eqnarray} \tag{ 75 }$$

Integration of this equation results in

$$\begin{eqnarray}&&{v}_{1}(t)=\displaystyle \frac{1}{({p}_{10}-1)| {v}_{00}| ({t}_{0}-t)+\tfrac{1}{{v}_{10}}}\end{eqnarray} \tag{ 76 }$$

with v₁(t₀) = v₁₀. For t → t_crit, the system develops a blow-up, where

$$\begin{eqnarray}&&{t}_{\mathrm{crit}}=\displaystyle \frac{1}{{v}_{10}({p}_{10}-1)| {v}_{00}| }+{t}_{0}.\end{eqnarray} \tag{ 77 }$$

To summarize, the choice of p₁(v₁) determines the temporal course of the flow and we have shown two contrasting examples, namely a decaying perturbation (for ${p}_{1}={\rm{const.}}$) and a perturbation that causes a blow-up of the system (for ${p}_{1}={p}_{10}{v}_{1}^{2}$). More complex choices of p₁(v₁), e.g., with higher orders of the polynomial, can lead to correspondingly diverse results such as multiple blow-ups and decays. Other time courses are quite possible and will occur, but one should keep in mind that the blow-up solutions can be damped by normal and anomalous dissipative processes, leading to a regular time-dependent behavior. However, the blow-up solutions indicate that nonlinear instabilities can develop in quasi-ideal systems and hence can lead to abrupt changes in the typical temporal scales.

The environments of certain types of evolved massive stars, such as the B[e] supergiants, display indications for disks, or multiple ring-like structures of yet unknown origin, which cannot be reconciled with the classical theory of a (viscous) outflowing disk (e.g., Lee et al. 1991; Okazaki 2001; Kurfürst et al. 2018). These configurations can be either steady (Kraus et al. 2016, 2023) or display temporal variability (Maravelias et al. 2018; Torres et al. 2018). Therefore, we have two problems, namely, (i) either long timescales (equilibria), or systems close to or converging to equilibrium states, or (ii) we need short timescales to explain the sudden appearance of new structures (such as detected by Oksala et al. 2012) resulting from nonlinear instabilities (occurring either within the circumstellar matter or already in the stellar atmosphere), which might be connected to collapse processes (finite-time singularities). To study the high diversity of geometrical structures and their possible formation mechanism we recourse to elementary mathematical tools of classical HD such as linear and nonlinear potential theory.

For the first case of long timescales, complex geometries of streamlines are known and can be constructed from the potential representation of equilibrium fields, e.g., potential fields (Laplace equation), or by nonlinear solutions of Grad–Shafranov-type equations, where closed and open streamline configurations (arcs, rings, radial structures) are possible (see, e.g., Nickeler et al. 2013, 2014). In our analysis, we include flows and gravity and find that the equations are of similar structure, which means that they can be solved in an analogous manner. We could prove that for a non-pressureless gas we can find quasi-Keplerian or even Keplerian rotation of the circumstellar matter.

To bridge the gap from these equilibrium structures to short timescales a homotopic modulation, i.e., a dilatation of the quasi-equilibria is performed, leading eventually to restrictions for the equilibrium values and to ordinary differential equations as constraints for the time dependencies of the fluid equations. These ordinary differential equations in time lead either to regular (decay or growth, or relaxation) or non-regular time dependencies where the latter includes finite-time singularities, which are known from (magneto-)hydrodynamical systems (e.g., Klapper et al. 1996; Nickeler & Karlický 2008), and a broad overview about mathematical solution techniques for different physical problems can be found in Galaktionov & Svirshchevskii (2007, and references therein).

For future investigations, non-separable time dependencies can be taken into account to get a better understanding of the variability and changes of the circumstellar matter. Moreover, for studies of the conditions under which rings and arcs can form within the circumstellar material it will be important to analyze in more detail the influence of different conformal mappings or multipole components, being able to create arcs, rings, and spiral-arm-like structures. Such an analysis will then allow us to derive the properties of the gas (density, temperature, and emissivity) which can be compared to observed quantities.

We wish to thank Prof. Dr. Wolfgang Glatzel for the warm hospitality during our visit to the Institut für Astrophysik (IAG), Georg-August-Universität Göttingen, and for the many fruitful discussions. This research made use of the NASA Astrophysics Data System (ADS). We acknowledge financial support from the Czech Science Foundation (GAČR, grant No. 20-00150S). The Astronomical Institute Ondřejov is supported by the project RVO:67985815. This project has received funding from the European Union's Framework Program for Research and Innovation Horizon 2020 (2014–2020) under the Marie Skłodowska-Curie grant agreement No. 823734.