††thanks: corresponding author

#  Entanglement Wedge Minimum Cross-Section for Holographic Aether Gravity

Chong-Ye Chen 1<EMAIL_ADDRESS>Wei Xiong 1
<EMAIL_ADDRESS>Chao Niu 1<EMAIL_ADDRESS>Cheng-Yong Zhang 1
<EMAIL_ADDRESS>Peng Liu 1<EMAIL_ADDRESS>1 Department of
Physics and Siyuan Laboratory, Jinan University, Guangzhou 510632, China

###### Abstract

We study the entanglement wedge cross-section (EWCS) in holographic Aether
gravity theory, a gravity theory with Lorentz symmetry violation while keeping
the general covariance intact. We find that only a limited parameter space is
allowed to obtain a black brane with positive Hawking temperature. Subject to
these allowed parameter regions, we find that the EWCS could exhibit non-
monotonic behaviors with system parameters. Meanwhile, the holographic
entanglement entropy (HEE), and the corresponding mutual information (MI), can
only exhibit monotonic behaviors. These phenomena suggest that the EWCS could
capture much more rich content of the entanglement than that of the HEE and
the MI. The role of the Lorentz violation in determining the behaviors of
quantum information-related quantities is also analyzed.

###### Contents

  1. I Introduction
  2. II Charged static solutions of Einstein-aether theory
    1. II.1 Introduction to holographic Aether Gravity
    2. II.2 Calculation of allowed parameter region
  3. III The Holographic entanglement entropy
  4. IV The Holographic mutual information
  5. V The Entanglement wedge minimum cross-section
  6. VI Discussion

##  I Introduction

Quantum information, as a unique feature of quantum systems, is gaining
increasing attention from areas such as holographic duality theory and
condensed matter theory. Many measures have been adopted to study quantum
entanglement, including entanglement entropy (EE), mutual information (MI),
Rényi entropy, and so on. According to their definitions, they can describe
different degrees of freedom of quantum systems. Among them, EE is the most
well-known physical quantity to capture the entanglement of a pure state.
However, EE is not suitable for describing the entanglement of mixed states.
Some other entanglement-related quantities, such as the MI, the Rényi
entanglement entropy, entanglement of purification (EoP), and reflected
entropy can be better candidate measures of mixed state entanglement
vidal:2002 ; Horodecki:2009review . However, one of the main problems with
quantum entanglement is that it is notoriously difficult to calculate.

Recently, gauge/gravity duality has been widely used to study strongly
correlated physics Ling:2015exa ; Landsteiner:2019kxb ; Ling:2014laa ;
Ling:2014bda . In addition to its power in strongly correlated systems, such
as condensed matter theory and QCD theory, gauge/gravity duality also
associates the information-related physical quantities in strongly correlated
systems with geometric quantities in dual gravitational systems Sakai:2004cn ;
Donos:2012js . The entanglement entropy of the dual quantum field theory has
been proposed proportional to the area of the minimum surface in the dual bulk
geometry, which was dubbed as the holographic entanglement entropy (HEE)
Ryu:2006bv . After that, many other holographic duals of quantum information-
related physical quantities have been proposed and studied Nishioka:2006gr ;
Klebanov:2007ws ; Pakman:2008ui ; Zhang:2016rcm ; Zeng:2016fsb ; Kudler-
Flam:2020url ; BabaeiVelni:2020wfl ; Sahraei:2021wqn ; Khoeini-
Moghaddam:2020ymm ; KumarBasak:2020eia . For example, the Rényi entropy has
been proposed proportional to the area of the minimum cosmic brane
Dong:2016fnf , which back-reacts on the background geometry. Meanwhile, the
entanglement of purification, reflected entropy and the entanglement
negativity have all been related to the minimum cross-section in the
entanglement wedge (EWCS) Takayanagi:2017knl ; Nguyen:2017yqw ; Kudler-
Flam:2018qjo ; Kusuki:2019zsp ; Dutta:2019gen ; Gong:2020pse ; Fu:2020oep ;
Liu:2021stu ; Zhang:2021edm ; Tamaoka:2018ned . In addition, quantum
complexity has been proposed proportional to the volume or the action of a
certain region, that can capture quantum information other than the
entanglement Susskind:2014rva . Furthermore, the butterfly velocity that
measures the speed of the propagation of quantum information has been found
related to the horizon of the black hole Shenker:2013pqa ; Sekino:2008he ;
Maldacena:2015waa ; Donos:2012js ; Blake:2016wvh ; Blake:2016sud ;
Ling:2016ibq ; Ling:2016wuy ; Wu:2017mdl ; Liu:2019npm . All of the above
progress become the cornerstone of the study of quantum information of the
strongly related systems in the framework of the holographic duality.

The Lorentz invariance is one of the principles of general relativity.
However, the violation of Lorentz invariance in the condensed matter systems
is common. To study the condensed matter systems with more practical
significance, the violation of Lorentz invariance in holographic duality
theory has been widely studied and found to exist in gravitational systems
such as massive gravity theories Vegh:2013sk ; Blake:2013bqa ; Blake:2013owa ;
Davison:2013jba ; deRham:2014zqa ; Baggioli:2014roa ; Alberte:2015isw .
However, a more natural approach is to explicitly break the Lorentz invariance
in the gravity theory. Recently, the theory of Aether gravity has been
proposed and studied as a kind of gravity theory that breaks the Lorentz
symmetry but preserves the general covariance Kai:2015nlbh ; Kai:2017ea .
Though HEE has been widely studied in many different holographic theories
Ling:2015dma ; Ling:2016wyr , the properties of mixed states entanglement such
as MI and EWCS in many holographic theories are still unclear. The main
purpose of this paper is to investigate the effect of Lorentz symmetry
violation on the entanglement of mixed states by studying EWCS, MI, and HEE in
the holographic Aether gravity theory.

In section II we will introduce the AdS Aether gravity model. Then we discuss
the properties of HEE (III), MI (IV) and EWCS (V). Lastly, we give a summary
in VI.

##  II Charged static solutions of Einstein-aether theory

###  II.1 Introduction to holographic Aether Gravity

The action of the $n$-dimensional Einstein-Aether-Maxwell theory reads
Kai:2017ea ,

$S=\int d^{n}x\frac{\sqrt{-g}}{16\pi
G_{ae}}(R-2\Lambda+\mathcal{L}_{ae}-\alpha F_{\mu\nu}F^{\mu\nu}),$ (1)

and

$\displaystyle\mathcal{L}_{ae}=$ $\displaystyle
c_{1}(\nabla_{\mu}u_{\nu})(\nabla^{\mu}u^{\nu})+c_{2}(\nabla^{\nu}u_{\nu})^{2}+c_{3}(\nabla_{\nu}u_{\mu})(\nabla^{\mu}u^{\nu})$
(2) $\displaystyle-
c_{4}u^{\mu}u^{\nu}(\nabla_{\mu}u_{\rho})(\nabla_{\nu}u^{\rho})+\lambda(u_{\nu}u^{\nu}+1).$

$F=dA$, and $A$ is the Maxwell field. $G_{ae}$ in (1) is the constant related
to Newton’s gravitational constant $G_{N}$ by $G_{ae}=(1-c_{14}/2)G_{N}$,
where $c_{ij}\equiv c_{i}+c_{j}$. The $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$
the coupling constants, and $\lambda$ the Lagrange multiplier such that the
aether vector $u_{a}$ satisfies the timelike constraint $u_{a}u^{a}=-1$. The
aether vector $u_{a}$ can also be expressed as the one-form of a scalar field
$\phi$,

$u_{a}\equiv\frac{-\partial_{a}\phi}{\sqrt{-g^{bc}\partial_{b}\phi\partial_{c}\phi}},$
(3)

which is called the khronon scalar Kai:2017ea ; Kai:2015nlbh . The speed of
the khronon scalar, defined as the speed of the mode of the perturbation of
$\phi$, is given by

$c_{\phi}^{2}=\frac{c_{123}}{c_{14}}.$ (4)

Eq. (4) shows that the most significant cases are $c_{123}=0$
$(c_{\phi}^{2}=0)$ and $c_{14}=0$ $(c_{\phi}^{2}\to\infty)$. Specifically, we
choose the case $c_{14}=0$ in this work. The ansatz of solution is given by

$\displaystyle ds^{2}$
$\displaystyle=-F(r)dv^{2}+2drdv+r^{2}h_{ij}dx^{i}dx^{j}$ (5) $\displaystyle
u^{a}$ $\displaystyle=u^{v}(r)(e_{v})^{a}+V(r)(e_{r})^{a}$ $\displaystyle
A_{a}$ $\displaystyle=A_{0}(r)(e_{v})_{a}$

where $u^{v}(r)=(V(r)+\sqrt{F(r)+V(r)^{2}})/F(r)$. The form of
$h_{ij}dx^{i}dx^{j}$ is given by

$h_{ij}dx^{i}dx^{j}=\begin{cases}d\theta^{2}+\sin^{2}\theta
d\Omega_{n-3}^{2}&k=1\\\ dx_{i}dx^{i}&k=0\\\ d\theta^{2}+\sinh^{2}\theta
d\Omega_{n-3}^{2}&k=-1\end{cases}\ i,j=1,2\ \cdots\ n-2.$ (6)

There are three options, sphere, planar or hyperbolic spacetime metric that
can be used as induced metric, corresponding to $k=1,0,-1$. Specifically, we
choose $k=0$, and induced metric component $h_{ij}$ takes the form of
$\delta_{ij}$. In the case of $c_{14}=0$, $n=4$ and $k=0$, we follow from
Kai:2017ea  and the solution is given by

$\displaystyle A_{0}(r)$ $\displaystyle=A_{c}-\frac{Q}{r}$ (7) $\displaystyle
F(r)$
$\displaystyle=\frac{2M_{a}}{r}+\frac{Q^{2}}{r^{2}}+\frac{4c_{13}B^{2}}{r^{4}}-\frac{2\Lambda-\rho_{a}B^{2}r_{z}^{-6}}{6}r^{2}$
$\displaystyle V(r)$ $\displaystyle=2Br(r_{z}^{-3}-r^{-3})$

where $\rho_{a}=24c_{13}-3c_{123}$. $A_{c}$, $B$ and $r_{z}$ are constants.
$M_{a}$ and $Q$ are the mass and charge of the black hole respectively.
$\Lambda$ is the cosmological constant. For convenience, we let
$r_{z}\to\infty$, $\beta\equiv c_{13}B^{2}$ and $\Lambda=-3$. The function $F$
in (7) becomes

$\displaystyle F(r)=$
$\displaystyle\frac{2M_{a}}{r}+\frac{{Q}^{2}}{r^{2}}+\frac{4\beta}{r^{4}}+r^{2}.$
(8)

For later calculation of the minimum surface, we transform the coordinate
system and the metric takes the form

$ds^{2}=-F(r)dt^{2}+F^{-1}(r)dr^{2}+r^{2}dx_{i}dx^{i},\quad i=1,2.$ (9)

At the horizon $r=r_{h}$ we have $F(r_{h})=0$. We can see that $r_{h}$
determines $M_{a}$ and the Hawking temperature reads,

$T=\frac{F^{\prime}(r_{h})}{4\pi}=-\frac{12\beta+Q^{2}r_{h}^{2}-3r_{h}^{6}}{4\pi
r_{h}^{5}},$ (10)

The system can be specified by $(T,\beta,Q)$. In addition, the aether vector
induces the Lorentz violation, and from the aether action (2) we see that
$c_{1},\,c_{2},\,c_{3},\,c_{4}$ associate with it. However, for the solution
that we are addressing, the overall effect of Lorentz violation is only
reflected in the parameter $\beta$, which we name as the Lorentz violation
parameter.

###  II.2 Calculation of allowed parameter region

A reasonable black hole system should have a non-negative Hawking temperature
$T$ and a non-negative horizon radius $r_{h}$. In addition, $Q$ always appears
as a squared term, so we only discuss the case $Q\geqslant 0$ here. In
summary, we discuss the physically reasonable parameter regions satisfying,

$T\geqslant 0,\;r_{h}\geqslant 0,\;Q\geqslant 0.$ (11)

The Hawking temperature (10) suggests that not any choice of $\beta$ satisfies
(11), hence we numerically work out the allowed regions. Firstly, for
$\beta<0$, the relationship between $T$ and $r_{h}$, $Q^{2}$ and $r_{h}$ are
two non-monotonic functions, and $T$ and $Q$ can have lower bounds (see Fig.
1).

Figure 1: The left plot: $T$ vs $r_{h}$ with fixed $Q$ and $\beta$.
Apparently, $T$ has a lower bound here. The right plot: $Q$ versus $r_{h}$
with fixed $T$ and $\beta$. $Q$ has a lower bound here.

Then we can get the minimum values of $Q$ or $T$ by locating the extreme
points. These extreme points are also called critical parameters, and the
surface formed by them is the boundary between the allowed parameter regions
and the non-physical parameter regions. For $\beta>0$, however, the situation
is different. It can be derived from the Hawking temperature (10) that if the
root $r_{h}$ is very large, the corresponding $T$ will also be very large;
while when the root $r_{h}$ tends to $0$, the corresponding $T$ will be
negative infinity. Next, we show the numerical results of the allowed
parameter regions.

We use NSolve in Mathematica to solve the critical values of $(Q,\beta,T)$,
then we obtain three contour plots (see Fig. 2).

Figure 2: From left to right are the contour plots of critical temperature,
critical $Q$ and critical $\beta$.

These plots show the critical value of $T,\,Q$, and $\beta$ as functions of
the other two parameters, respectively. More specifically, the critical $T$
increases with decreasing $Q$ and $\beta$; the critical $Q$ increases with
decreasing $\beta$ and $T$; the critical $\beta$ increases with decreasing $Q$
and $T$. In addition, the red curves in the left and middle plot of Fig. 2 is
the boundary where the critical $T$ or $Q$ vanishes. In general, the critical
parameter of each in $(Q,\beta,T)$ decreases along the direction of the other
two parameters. Moreover, as mentioned above, when $\beta>0$, the critical
value of both $Q$ and $T$ will drop to $0$. Bounded by these regions, we
calculate and discuss the properties of HEE, MI and EWCS in the section III,
IV and V.

##  III The Holographic entanglement entropy

One of the most important features distinguishing quantum systems from
classical ones is entanglement, which can be measured by many physical
quantities. As the most well-known measure, EE depicts the entanglement
between a subsystem and its complement. Given a system composed of disjoint
$A$ and $B$, the subsystem $A$ is described by a reduced density matrix
$\rho_{A}=\text{Tr}_{B}\rho_{\text{total}}$. To characterize the entanglement
between $A$ and $B$, EE is defined as the von Newmann entropy of the reduced
density matrix,

$S_{A}(|\psi\rangle)=-\text{Tr}\left[\rho_{A}\log\rho_{A}\right],\quad\rho_{A}=\text{Tr}_{B}\left(|\psi\rangle\langle\psi|\right).$
(12)

This definition immediately leads to $S_{A}=S_{B}$ for pure states
Chuang:2002book . Though EE has been widely considered a good entanglement
measure for pure states, it is not suitable for describing the entanglement of
mixed states. Because even if the degrees of freedom of $A$ and $B$ are not
entangled, such as the direct product states, they can still have non-zero EE.
Many new entanglement measures have been proposed to characterize the
entanglement of mixed states, among which MI is the most commonly used one
vidal:2002 ; Horodecki:2009review . In holographic duality theory, the EE was
associated with the area of the minimum surface stretching into the bulk of
the dual gravity systems (see the left plot of Fig. 3) Ryu:2006bv .

$x$$y$$z$

$x$$y$$z$

Figure 3: The left plot: The minimum surface for a given width $w$. The right
plot: The minimum cross-section (green surface) of the entanglement wedge.

In this paper, we focus on the partition of infinite strips along
$y$-direction, where solving the minimum surfaces only involves ordinary
differential equations. We parametrize the minimum surface with the angle
$\theta\equiv\arctan(z/x)$ (see Fig. 10), which can facilitate the solving of
the minimum surface Liu:2020blk ; Liu:2021rks . The first step of this
numerical method is to discretize the angle with Gauss-Lobatoo collocation
Boyd:2001 . Because of the nonlinearity of the equation of motions for the
minimum surfaces, we must apply the Newton-Raphson iteration method to find
out the minimal surface. Equipped with this numerical method, we study the
properties of HEE in the Aether gravity model.

First, we show the HEE versus $T$ in Fig. 4, from which we can find that HEE
increases with $T$, whether the $\beta$ is negative or positive.

Figure 4: HEE vs $T$ in different value of the $\beta$ when $Q=4$ and width
$w=1$.

The temperature behavior of the HEE depends on the relation between the
$r_{h}$ and the temperature. When the temperature increases, the $r_{h}$
increases, hence the minimum surface tends to approach the horizon. As a
consequence, the HEE increases with the temperature. Physically speaking, this
is as expected since the entanglement entropy becomes more and more dominated
by the thermal entropy, which monotonically increases with the temperature.

Next, we show the HEE versus $Q$ in Fig. 5, from which we can find that HEE
increases with $Q$ and width $w$.

Figure 5: HEE vs $Q$ in different value of width $w$ when $\beta=1$ and
$T=0.3$.

Similar to the temperature behavior, the HEE behavior with $Q$ can also be
attributed to the increase of $r_{h}$ when increasing $Q$. Similar behaviors
have also been obtained in AdS-RN system Liu:2019qje .

The particularly interesting phenomenon is the HEE along the Lorentz violation
parameter $\beta$ (Fig. 6).

Figure 6: HEE vs $\beta$ in different value of width $w$ when $Q=3$ and
$T=0.2$. The HEE increases with the increase of $\beta$.

The underlying reason for the increasing HEE with $\beta$ is also the same as
that of the $T$ and $Q$. The partial derivative of $r_{h}$ with $\beta$ can be
obtained as

$\partial_{\beta}r_{h}=\frac{4r_{h}}{20\beta+Q^{2}r_{h}^{2}+r_{h}^{6}},$ (13)

which shows that $\partial_{\beta}r_{h}>0$ for physically allowed regions. The
horizon radius of the black brane becomes larger as the Lorentz violation
effect becomes more significant. Since the Lorentz violation effect is induced
by the introduction of time-like field $u^{a}$, it can be clearly seen from
Eq. (2) and Eq. (7) that the increase of horizon radius of black brane is
caused by the increasing of the action of the matter field $u^{a}$.
Furthermore, it leads to the increase of the entanglement entropy of the dual
boundary theory.

Although they have the same monotonicity, it is worth noting that there are
obvious differences in the details of HEE behaviors with $Q$, $T$, and
$\beta$. First, compared with the increasing behavior of HEE with $\beta$, HEE
increases more rapidly with $T$ and $Q$. This phenomenon can be understood
from Eq. (13), where $\partial_{\beta}r_{h}$ decreases with the increase of
$\beta$. Therefore, with the increase of $\beta$, HEE will increase, meanwhile
the slope decreases.

After elaborating on the properties of HEE, we now discuss the properties of
MI, the entanglement measure of mixed states, in the Aether gravity.

##  IV The Holographic mutual information

The MI for disjoint $A\cup B$ is defined as

$I\left(A,B\right):=S\left(A\right)+S\left(B\right)-S\left(A\cup B\right),$
(14)

When $\rho_{AB}=\rho_{A}\otimes\rho_{B}$, it can be verified that
$I\left(A,B\right)=0$. Therefore, MI can recognize that the product states are
not entangled. Based on the numerical calculation of HEE in the previous
section, we can calculate the MI in aether gravity, as the definition of MI is
directly related to HEE.

Figure 7: The MI vs $T$ in different value of $T$. It shows the change rule
between the MI and two parameters $T$ and $\beta$ that the MI decreases with
the increase of $T$ and $\beta$.  Figure 8: The MI vs $\beta$ in different
configurations $(a,b,c)$ by choosing $a=0.6$, $b=0.1$ with different value of
$c$.The MI decreses with the increasing $\beta$. However, MI increases with
the increasing of c.This means that MI increases when the configuration tends
to be symmetrical. Figure 9: The MI vs $Q$ in different configurations
$(a,b,c)$ taking the same proportion between $a$ ,$b$ and $c$. The MI
decreases with the increasing charge $Q$ and the larger $(a,b,c)$.

We demonstrate the behavior of MI in the Aether gravity as follows. In
general, given a configuration $(a,b,c)$, the behavior of MI is opposite to
that of HEE - MI decreases with the increase of $T$, $Q$ and $\beta$. First,
we can see from Fig. 7, Fig. 8 and Fig. 9 that MI always decreases with the
increase of $T,\,\beta,\,Q$. The reason responsible for the opposite monotonic
behavior of the MI to that of the HEE directly results from the HEE. According
to the definition (14), a non-trivial MI will be $I=S_{A}+S_{C}-S_{B}-S_{A\cup
B\cup C}$, the last term $S_{A\cup B\cup C}$ is more affected by the
deformation caused by the deviation from the AdS, and therefore the MI shows
the opposite behavior.

Moreover, we can find from Fig. 8 that MI increases with the increase of $c$
when fixing $a$ and $b$. The reason behind this phenomenon is that increasing
the size of the subregion can increase the degrees of freedom that can
entangle with other degrees of freedom. In addition, we discuss the behavior
of MI when varying the sizes of $a$, $b$, and $c$ in equal proportions (see
Fig. 9), where different curve corresponds to different value of $(a,b,c)$.
The MI decreases when uniformly increasing the values of $(a,b,c)$. This shows
that the decrease of entanglement caused by the increase of the separation
plays a major role when increasing the subregions and the separation in the
same proportion. In fact, this is a reasonable result. Because the
entanglement between the degrees of freedom usually decays more rapidly with
the increase of the separation, when the separation increases to a certain
value, the subregions will be disentangled.

The above phenomena show that MI is directly determined by HEE, thus MI may
not be a good measure of mixed state entanglement. Therefore we need to resort
to other mixed state entanglement measures. In the next section, we will study
a new entanglement measure, the EWCS, in the Aether gravity.

##  V The Entanglement wedge minimum cross-section

The entanglement wedge minimum cross-section (EWCS) has been associated with
several different mixed state entanglement measures, such as entanglement of
purification, the reflected entropy, odd entropy, and so on. It involves the
purification process of mixed states. When the entanglement wedge exists
(i.e., where MI is non-trivial), the mixed entanglement measures
$E_{W}\left(\rho_{AB}\right)$ is proportional to the area of the minimum
cross-section $\Sigma_{AB}$ in it Takayanagi:2017knl ,

$E_{W}\left(\rho_{AB}\right)=\min_{\Sigma_{AB}}\left(\frac{\text{Area}\left(\Sigma_{AB}\right)}{4G_{N}}\right).$
(15)

EWCS vanishes when the entanglement wedge becomes disconnected, i.e., when MI
vanishes.

Though the explicit definition is given, the following three reasons explain
the difficulty of solving the EWCS. First, the equations of motion of minimum
surfaces are highly nonlinear, which are usually hard to solve. Second, the
minimum cross-section is hard to find because each cross-section itself is
already a local minimum, which means that minimizing the cross-section is a
second-order minimization. Last but not least, the coordinate singularity at
the horizon will easily sabotage the numerical precision, adding to the
difficulty of our numerics.

An efficient algorithm for solving EWCS has been developed by using the
boundary condition that the minimum cross-section must be locally orthogonal
to the entanglement wedge Liu:2020blk . Fig. 10 shows a schematic diagram of
the main ideas of the algorithm for solving EWCS.

Figure 10:  Schematic diagram of solving EWCS. The inner blue curve
$C_{1}(\theta_{1})$ and the outer orange curve $C_{2}(\theta_{2})$ in the
figure are the corresponding minimum surfaces with widths of $b$ and $a+b+c$
respectively. The red curve is the minimum surface connecting $p_{1}$ on
$C_{1}$ and $p_{2}$ on $C_{2}$. There are blue and purple arrows at $p_{1}$
and $p_{2}$, representing the tangent vector along the red curve and the
tangent vector along $C_{1}$ or $C_{2}$. The horizontal black dashed line
represents the horizon of the black brane.

We focus on the EWCS of parallel infinite strips in a homogeneous background,

$ds^{2}={g_{tt}}dt^{2}+g_{zz}dz^{2}+g_{xx}dx^{2}+g_{yy}dy^{2},$ (16)

where the asymptotic AdS boundary resides at $z=0$. The homogeneity implies
that $z$ is the only variable of the metric components $g_{\mu\nu}$. For a
biparty subsystem with minimum surfaces
$C_{1}(\theta_{1}),\,C_{2}(\theta_{2})$, we work out the minimum surface
$C_{p_{1},p_{2}}$ connecting $p_{1}\in C_{1}$ and $p_{2}\in C_{2}$. By
parametrizing $C_{p_{1},p_{2}}$ with $z$, the area of $C_{p_{1},p_{2}}$ reads,

$A=\int_{C_{p_{1},p_{2}}}\sqrt{g_{xx}g_{yy}x^{\prime}(z)^{2}+g_{zz}g_{yy}}dz.$
(17)

The equation of motion from minimizing $A$ reads,

$x^{\prime}(z)^{3}\left(\frac{g_{xx}g_{yy}^{\prime}}{2g_{yy}g_{zz}}+\frac{g_{xx}^{\prime}}{2g_{zz}}\right)+x^{\prime}(z)\left(\frac{g_{xx}^{\prime}}{g_{xx}}+\frac{g_{yy}^{\prime}}{2g_{yy}}-\frac{g_{zz}^{\prime}}{2g_{zz}}\right)+x^{\prime\prime}(z)=0,$
(18)

with boundary conditions,

$x(z(\theta_{i}))=x(\theta_{i}),\quad i=1,2.$ (19)

The local orthogonal relation between the minimum cross-section and the
entanglement wedge leads to,

$\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial\theta_{1}}\right\rangle_{p_{1}}=0,\quad\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial\theta_{2}}\right\rangle_{p_{2}}=0,$ (20)

where $\langle\cdot,\cdot\rangle$ represents the vector product under the
metric $g_{ab}$. In order to control the numerical precision, we adopt the
normalized local orthogonal relation,

$Q_{1}(\theta_{1},\theta_{2})\equiv\left.\frac{\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial\theta_{1}}\right\rangle}{\sqrt{\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial
z}\right\rangle\left\langle\frac{\partial}{\partial\theta_{1}},\frac{\partial}{\partial\theta_{1}}\right\rangle}}\right|_{p_{1}}=0,\quad
Q_{2}(\theta_{1},\theta_{2})\equiv\left.\frac{\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial\theta_{1}}\right\rangle}{\sqrt{\left\langle\frac{\partial}{\partial
z},\frac{\partial}{\partial
z}\right\rangle\left\langle\frac{\partial}{\partial\theta_{2}},\frac{\partial}{\partial\theta_{2}}\right\rangle}}\right|_{p_{2}}=0.$
(21)

Now, the EWCS can be given when we find the minimum surface on which
$(\theta_{1},\theta_{2})$ of the endpoints satisfy (21). To this end, we
variate the endpoints satisfying the local perpendicular conditions by the
Newton-Raphson method. Next, we study the relationship between the Aether
gravity and the EWCS based on the above techniques.

The EWCS behaves distinctly from the HEE or MI in Aether gravity. We fix the
configuration $(a,b,c)$ to $(0.6,0.1,0.3)$ since the main phenomenon is
independent of this. First, the EWCS shows non-monotonicity with the increase
of $Q$ (see Fig. 11).

Figure 11: The EWCS vs $Q$ in different value of $\beta$ and the configuration
$(a,b,c)$ is $(0.6,0.1,0.3)$. The EWCS behave non-monotonicity with the
increase of charge $Q$, which also will vanish if $\beta$ takes a rather
larger value. Three red points represent the critical charge $Q$.

EWCS decreases with the increase of $Q$ at the beginning and then starts to
increase. However, when $\beta$ increases to $0.15$, we can see that EWCS
increases monotonically with the increase of $Q$. The non-monotonicity of it
will vanish if the $\beta$ takes a larger value. In addition, three red points
$p_{1}$, $p_{2}$ and $p_{3}$ in Fig. 12 show the critical charge $Q$ of the
cases $\beta=-0.15,-0.05,-0.025$ in Fig. 11 subject to the allowed parameter
region shown in Fig. 2.

Figure 12: Critical charge $Q$ in case of $T=0.25$

For EWCS behavior along $\beta$-direction, the EWCS first decreases and then
increases with the increasing $\beta$ (see Fig. 13).

Figure 13: The EWCS vs $\beta$ in different value of $Q$ and the configuration
$(a,b,c)$ is $(0.6,0.1,0.3)$.The EWCS shows the non-monotonicity with the
increase of $\beta$ at every $Q$. We can see that the EWCS first decreases
with $Q$ and $\beta$, and then turns into an increase.

Also, the intersections shown in Fig. 13 again reflects the non-monotonicity
of EWCS with $Q$. When $\beta$ is small, EWCS decreases rapidly with the
growth of $\beta$, which can also be understood through the relationship
between $r_{h}$ and $\beta$. It is noted that $r_{h}$ increases rapidly with
the increase of $\beta$ when $\beta$ is small. When $r_{h}$ increases, the
minimum surface will be close to the horizon of the black hole and moves away
from the AdS boundary. Because the area of EWCS is mainly contributed by AdS
boundary, thus when $\beta$ is small, EWCS will decrease rapidly with the
increase of $\beta$. With the increase of $\beta$, $r_{h}^{\prime}(\beta)$
becomes smaller, so the growth of EWCS slows down gradually. Therefore, EWCS
will present a relatively flat area. When $\beta$ is large, the increasing
behavior of EWCS with the increase of $\beta$ is no longer controlled by the
deviation from AdS. Instead, we must consider the specific contribution of
bulk geometry. Specifically, the variation of the EWCS comes from the
variation of the background metric and the minimum surface,

$\displaystyle\delta E_{W}$
$\displaystyle=\delta\int_{C_{p_{1},p_{2}}}\sqrt{g_{yy}\left(g_{xx}dx^{2}+g_{rr}dr^{2}\right)}$
(22) $\displaystyle=\int_{C_{p_{1},p_{2}}}\frac{\delta E_{W}}{\delta
g_{\mu\nu}}\delta g_{\mu\nu}+\int_{C_{p_{1},p_{2}}}\frac{\delta E_{W}}{\delta
C_{p_{1},p_{2}}}\delta C_{p_{1},p_{2}}.$

The $\frac{\delta E_{W}}{\delta C_{p_{1},p_{2}}}$ of (22) is the equation of
motion, which must vanish. The first term on the right-hand side of (22)
denotes the contribution to the EWCS from the variation of the background
metric, which results from the variation of the $F_{a}(r)$.

$\displaystyle\partial_{\beta}E_{W}$
$\displaystyle=\int_{C_{p_{1},p_{2}}}\partial_{\beta}g_{rr}\sqrt{\frac{g_{yy}}{g_{xx}dx^{2}+g_{rr}dr^{2}}}dr^{2}$
(23)
$\displaystyle=\int_{C_{p_{1},p_{2}}}\left(-\frac{\partial_{\beta}F_{a}(r)}{F_{a}^{2}(r)}\right)\sqrt{\frac{g_{yy}}{g_{xx}dx^{2}+g_{rr}dr^{2}}}dr^{2}.$

After replacing $M_{a}$ and $Q$ with horizon boundary conditions and Hawking
temperature, we obtain that,

$F_{a}(r)=-\frac{-r^{4}r_{h}+4\pi
Tr_{h}^{4}-3r_{h}^{5}+rr_{h}^{4}+Q^{2}r}{r^{2}r_{h}}-\frac{4\beta\left(3r^{2}r_{h}-r_{h}^{3}+r^{3}\right)}{r^{4}r_{h}^{3}}.$
(24)

Together with the expression $\partial_{\beta}r_{h}$ (13), we will find that,

$\partial_{\beta}F_{a}(r)=\frac{4}{5r^{4}r_{h}^{3}}\left(\frac{2r^{2}r_{h}^{2}\left(3r_{h}^{4}\left(-11r_{h}+2r+10\pi
T\right)-Q\left(3r_{h}+4r\right)\right)}{20\beta+Qr_{h}^{2}+r_{h}^{6}}-9r^{2}r_{h}+5r_{h}^{3}-2r^{3}\right).$
(25)

When $\beta$ is relatively large,

$\partial_{\beta}F_{a}(r)\simeq-\frac{4\left(9r^{2}r_{h}-5r_{h}^{3}+2r^{3}\right)}{5r^{4}r_{h}^{3}}<0,$
(26)

and hence the second term in (24) is negative and $\partial_{\beta}E_{W}>0$ in
(23). This arguments are also supported by numerics shown in Fig. 14, where
varying only the metric will render a increasing behavior of the EWCS with
$\beta$.

Figure 14:  The EWCS vs $\beta$ in different values of $Q$ and the
configuration $(a,b,c)$ is $(0.6,0.1,0.3)$ when fixing the minimum surface and
varying the background geometry.

From the dual picture, when the Lorentz symmetry violation effect is weak, the
mixed state entanglement decreases rapidly with the enhancement of Lorentz
symmetry violation; when the Lorentz symmetry violation effect is strong, the
mixed state entanglement of the system increases slowly when further enhancing
the Lorentz symmetry violation.

In Fig. 15, EWCS decreases with the increase of temperature. This is in line
with physical expectations. Usually, heating up a quantum system will destroy
the quantum entanglement between the two subregions. Therefore, EWCS behavior
with the temperature here will reflect this. Also, in the left plot of Fig. 15
the curves with different values of $\beta$ cross each other.

Figure 15: The former plot is EWCS vs $T$ in different values of $\beta$. The
EWCS decrease with the increasing temperature. There are several red points
and purple points that show the non-monotonicity of the EWCS vs $\beta$. The
latter plot is EWCS vs $T$ in different values of $Q$. The configuration
$(a,b,c)$ of them are both $(0.6,0.1,0.3)$

Specifically, we can see the red points and the purple points on curves, which
correspond to the points in Fig. 16, from which it is easy to see that EWCS
behaves non-monotonically with $\beta$ at a fixed temperature.

Figure 16: The EWCS vs $\beta$ in different value of $T$ in cases of $T=0.05$
and $T=0.2$.

The point $A_{1}$ is lower than the point $A_{5}$ in case of $T=0.5$, however
when changing the temperature to $T=0.2$, the point $B_{1}$ become higher than
point $B_{5}$. There are also several crossings of $E_{W}$ vs $T$ in the right
plot of Fig. 15, where different curves represent different values of $Q$. It
is also easy to see the non-monotonic relationship between the $E_{W}$ and the
$Q$.

In summary, the EWCS shows non-monotonicity along with charge $Q$ and the
Lorentz violation parameter $\beta$ in Aether gravity. Specifically, the EWCS
first decreases with the increase of them then increases with the increase of
them. However, EWCS decreases monotonically with temperature $T$.

##  VI Discussion

In this paper, we studied the properties of mixed state entanglement and
entanglement entropy in the Aether gravity theory with Lorentz symmetry
violation. HEE and MI are found to monotonically change with the charge $Q$,
the Lorentz violation parameter $\beta$, and temperature $T$. First, HEE
increases with the increasing $Q$, $\beta$, $T$, and the width $w$. Moreover,
as a measure of mixed state entanglement, MI shows exactly the opposite
monotonicity to HEE. These results are independent of the specific
configuration. More importantly, we found that EWCS behaves very differently
from HEE or MI in Aether gravity. With the increasing $T$, EWCS decreases
monotonically. However, EWCS shows non-monotonicity with $Q$ and $\beta$: it
first decreases with the increasing $Q$ and $\beta$ and then increases with
them. Based on the analytical treatments and numerical results, we show that
EWCS behaves non-monotonically in the direction of $\beta$ due to the special
role of Lorentz violation parameter $\beta$ in geometry. When the Lorentz
violation parameter of the system is small, EWCS decreases rapidly with the
increase of $\beta$. However, when the Lorentz violation parameter is large,
that is, the $\beta$ is very large, it has little influence on the
entanglement properties when further increasing $\beta$. These phenomena will
lay a foundation for the further study of the gravity system with Lorentz
violation.

One of the topics worthy of further study is to examine other gravity models
with Lorentz violation to check whether the effect of Lorentz violation is
consistent with those discussed in this paper. In addition, it is worth
exploring other models that do not have explicit Lorentz violation but have
obvious Lorentz violation in their dual systems. It is desirable to reveal the
properties of mixed state entanglement in these models to further understand
the influence of Lorentz violation on quantum information-related physical
quantities. The above research will help to understand the relationship
between systems with and without explicit Lorentz violation, and to further
understand the duality of Lorentz violation in condensed matter theory.

## Acknowledgments

Peng Liu would like to thank Yun-Ha Zha for her kind encouragement during this
work. This work is supported by the Natural Science Foundation of China under
Grant No. 11805083, 11905083, 12005077 and Guangdong Basic and Applied Basic
Research Foundation (2021A1515012374)

## References

  * (1) Vidal, G. and Werner, R.F., 2002. “A computable measure of entanglement”, Physical Review A, 65(3), p.032314. quant-ph:0102117
  * (2) Horodecki, R., Horodecki, P., Horodecki, M., Horodecki, K. (2009). “Quantum entanglement.” Reviews of modern physics, 81(2), 865.
  * (3) Y. Ling, P. Liu and J. P. Wu, “A novel insulator by holographic Q-lattices,” JHEP 1602, 075 (2016)
  * (4) K. Landsteiner, Y. Liu and Y. W. Sun, “Holographic topological semimetals,” Sci. China Phys. Mech. Astron. 63 (2020) no.5, 250001 [arXiv:1911.07978 [hep-th]].
  * (5) Y. Ling, P. Liu, C. Niu, J. P. Wu and Z. Y. Xian, “Holographic Superconductor on Q-lattice,” JHEP 1502, 059 (2015)
  * (6) Y. Ling, P. Liu, C. Niu, J. P. Wu and Z. Y. Xian, “Holographic fermionic system with dipole coupling on Q-lattice,” JHEP 1412, 149 (2014)
  * (7) T. Sakai and S. Sugimoto, Prog. Theor. Phys. 113, 843-882 (2005) doi:10.1143/PTP.113.843 [arXiv:hep-th/0412141 [hep-th]].
  * (8) A. Donos and S. A. Hartnoll, “Metal-insulator transition in holography”, Nature Phys. 9, 649 (2013) [arXiv:1212.2998].
  * (9) S. Ryu and T. Takayanagi, “Holographic derivation of entanglement entropy from AdS/CFT,” Phys. Rev. Lett. 96, 181602 (2006) [hep-th/0603001].
  * (10) T. Nishioka and T. Takayanagi, “AdS Bubbles, Entropy and Closed String Tachyons,” JHEP 0701, 090 (2007) [hep-th/0611035].
  * (11) I. R. Klebanov, D. Kutasov and A. Murugan, “Entanglement as a probe of confinement,” Nucl. Phys. B 796, 274 (2008) [arXiv:0709.2140 [hep-th]].
  * (12) A. Pakman and A. Parnachev, “Topological Entanglement Entropy and Holography,” JHEP 0807, 097 (2008) [arXiv:0805.1891 [hep-th]].
  * (13) S. J. Zhang, “Holographic entanglement entropy close to crossover/phase transition in strongly coupled systems,” Nucl. Phys. B 916, 304 (2017) [arXiv:1608.03072 [hep-th]].
  * (14) X. X. Zeng and L. F. Li, “Holographic Phase Transition Probed by Nonlocal Observables,” Adv. High Energy Phys. 2016, 6153435 (2016) [arXiv:1609.06535 [hep-th]].
  * (15) J. Kudler-Flam, Y. Kusuki and S. Ryu, “Correlation measures and the entanglement wedge cross-section after quantum quenches in two-dimensional conformal field theories,” JHEP 04 (2020), 074 [arXiv:2001.05501 [hep-th]].
  * (16) K. Babaei Velni, M. R. Mohammadi Mozaffar and M. H. Vahidinia, JHEP 08 (2020), 129 doi:10.1007/JHEP08(2020)129 [arXiv:2005.05673 [hep-th]].
  * (17) M. Sahraei, M. J. Vasli, M. R. M. Mozaffar and K. B. Velni, JHEP 08 (2021), 038 doi:10.1007/JHEP08(2021)038 [arXiv:2105.12476 [hep-th]].
  * (18) S. Khoeini-Moghaddam, F. Omidi and C. Paul, JHEP 02 (2021), 121 doi:10.1007/JHEP02(2021)121 [arXiv:2011.00305 [hep-th]].
  * (19) J. Kumar Basak, V. Malvimat, H. Parihar, B. Paul and G. Sengupta, [arXiv:2002.10272 [hep-th]].
  * (20) X. Dong, “The Gravity Dual of Rényi Entropy,” Nature Commun. 7, 12472 (2016) [arXiv:1601.06788 [hep-th]].
  * (21) T. Takayanagi and K. Umemoto, “Holographic Entanglement of Purification,” arXiv:1708.09393 [hep-th].
  * (22) P. Nguyen, T. Devakul, M. G. Halbasch, M. P. Zaletel and B. Swingle, “Entanglement of purification: from spin chains to holography,” JHEP 1801, 098 (2018) [arXiv:1709.07424 [hep-th]].
  * (23) J. Kudler-Flam and S. Ryu, “Entanglement negativity and minimal entanglement wedge cross-sections in holographic theories,” Phys. Rev. D 99 (2019) no.10, 106014 [arXiv:1808.00446 [hep-th]].
  * (24) Y. Kusuki, J. Kudler-Flam and S. Ryu, “Derivation of Holographic Negativity in AdS3/CFT2,” Phys. Rev. Lett. 123 (2019) no.13, 131603 [arXiv:1907.07824 [hep-th]].
  * (25) S. Dutta and T. Faulkner, “A canonical purification for the entanglement wedge cross-section,” JHEP 03 (2021), 178 [arXiv:1905.00577 [hep-th]].
  * (26) H. Gong, P. Liu, G. Fu, X. M. Kuang and J. P. Wu, “Informational properties of holographic Lifshitz field theory,” Chin. Phys. C 45, no.6, 065101 (2021) [arXiv:2009.00450 [hep-th]].
  * (27) G. Fu, P. Liu, H. Gong, X. M. Kuang and J. P. Wu, “Informational properties for Einstein-Maxwell-Dilaton Gravity,” [arXiv:2007.06001 [hep-th]].
  * (28) P. Liu and J. P. Wu, “Dynamic Properties of Two-Dimensional Latticed Holographic System,” [arXiv:2104.04189 [hep-th]].
  * (29) C. Y. Zhang, P. Liu, Y. Liu, C. Niu and B. Wang, “Evolution of Anti-de Sitter black holes in Einstein-Maxwell-dilaton theory,” [arXiv:2104.07281 [gr-qc]].
  * (30) K. Tamaoka, “Entanglement Wedge Cross Section from the Dual Density Matrix,” Phys. Rev. Lett. 122 (2019) no.14, 141601 [arXiv:1809.09109 [hep-th]].
  * (31) L. Susskind, “Computational Complexity and Black Hole Horizons,” [Fortsch. Phys. 64, 24 (2016)] Addendum: Fortsch. Phys. 64, 44 (2016)
  * (32) S. H. Shenker and D. Stanford, “Black holes and the butterfly effect,” JHEP 1403, 067 (2014) [arXiv:1306.0622 [hep-th]].
  * (33) Y. Sekino and L. Susskind, “Fast Scramblers,” JHEP 0810, 065 (2008) [arXiv:0808.2096 [hep-th]].
  * (34) J. Maldacena, S. H. Shenker and D. Stanford, “A bound on chaos,” JHEP 1608, 106 (2016) [arXiv:1503.01409 [hep-th]].
  * (35) M. Blake, “Universal Charge Diffusion and the Butterfly Effect in Holographic Theories,” Phys. Rev. Lett. 117, no. 9, 091601 (2016) [arXiv:1603.08510 [hep-th]].
  * (36) M. Blake, “Universal Diffusion in Incoherent Black Holes,” Phys. Rev. D 94, no. 8, 086014 (2016) [arXiv:1604.01754 [hep-th]].
  * (37) Y. Ling, P. Liu and J. P. Wu, “Holographic Butterfly Effect at Quantum Critical Points,” JHEP 1710, 025 (2017) [arXiv:1610.02669 [hep-th]].
  * (38) Y. Ling, P. Liu and J. P. Wu, “Note on the butterfly effect in holographic superconductor models,” Phys. Lett. B 768, 288 (2017) [arXiv:1610.07146 [hep-th]].
  * (39) S. F. Wu, B. Wang, X. H. Ge and Y. Tian, “Collective diffusion and quantum chaos in holography,” Phys. Rev. D 97, no. 10, 106018 (2018) [arXiv:1702.08803 [hep-th]].
  * (40) P. Liu, C. Niu and J. P. Wu, “The Effect of Anisotropy on Holographic Entanglement Entropy and Mutual Information,” Phys. Lett. B 796, 155 (2019) [arXiv:1905.06808 [hep-th]].
  * (41) D. Vegh, “Holography without translational symmetry,” [arXiv:1301.0537 [hep-th]].
  * (42) M. Blake and D. Tong, “Universal Resistivity from Holographic Massive Gravity,” Phys. Rev. D 88 (2013) no.10, 106004 [arXiv:1308.4970 [hep-th]].
  * (43) M. Blake, D. Tong and D. Vegh, “Holographic Lattices Give the Graviton an Effective Mass,” Phys. Rev. Lett. 112 (2014) no.7, 071602 [arXiv:1310.3832 [hep-th]].
  * (44) R. A. Davison, “Momentum relaxation in holographic massive gravity,” Phys. Rev. D 88 (2013), 086003 [arXiv:1306.5792 [hep-th]].
  * (45) C. de Rham, “Massive Gravity,” Living Rev. Rel. 17 (2014), 7 [arXiv:1401.4173 [hep-th]].
  * (46) M. Baggioli and O. Pujolas, “Electron-Phonon Interactions, Metal-Insulator Transitions, and Holographic Massive Gravity,” Phys. Rev. Lett. 114 (2015) no.25, 251602 [arXiv:1411.1003 [hep-th]].
  * (47) L. Alberte, M. Baggioli, A. Khmelnitsky and O. Pujolas, “Solid Holography and Massive Gravity,” JHEP 02 (2016), 114 [arXiv:1510.09089 [hep-th]].
  * (48) K. Lin, O. Goldoni, M. F. da Silva and A. Z. Wang “New look at black holes: Existence of universal horizons” [arXiv:1410.6678v3 [gr-qc]]
  * (49) K. Lin, F. H. Ho and W. L. Qian “Charged Einstein-æther black holes in n-dimensional spacetime” [arXiv:1704.06728v3 [gr-qc]]
  * (50) Y. Ling, P. Liu, C. Niu, J. P. Wu and Z. Y. Xian, “Holographic Entanglement Entropy Close to Quantum Phase Transitions,” JHEP 1604, 114 (2016)
  * (51) Y. Ling, P. Liu and J. P. Wu, “Characterization of Quantum Phase Transition using Holographic Entanglement Entropy,” Phys. Rev. D 93, no.12, 126004 (2016) [arXiv:1604.04857 [hep-th]].
  * (52) Nielsen, Michael A., and Isaac Chuang. “Quantum computation and quantum information.” (2002): 558-559.
  * (53) P. Liu and J. P. Wu, “Mixed state entanglement and thermal phase transitions,” Phys. Rev. D 104 (2021) no.4, 046017 [arXiv:2009.01529 [hep-th]].
  * (54) P. Liu, C. Niu, Z. J. Shi and C. Y. Zhang, “Entanglement Wedge Minimum Cross-section in Holographic Massive Gravity Theory,” [arXiv:2104.08070 [hep-th]].
  * (55) John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
  * (56) P. Liu, Y. Ling, C. Niu and J. P. Wu, “Entanglement of Purification in Holographic Systems,” JHEP 09 (2019), 071 [arXiv:1902.02243 [hep-th]].


# Accretion in the Binary System GG Carinae and Implications for B[e]
Supergiants

Amit Kashi1,2
1Department of Physics, Ariel University, Ariel, 4070000, Israel
2Astrophysics, Geophysics and Space Science (AGASS) Center, Ariel University,
Ariel, 4070000, Israel
E-mail<EMAIL_ADDRESS>
(Accepted 9 June 2023. Received 8 June 2023; In original form 19 April 2023)

###### Abstract

We simulate the hydrodynamics of the wind flow in the B[e] supergiant binary
system GG Carinae and obtain the mass accretion rate onto the secondary and
the observed lightcurve. We find an inhomogeneous Bondi-Hoyle-Lyttleton
accretion into a curved accretion tail, and confirm that the accretion rate is
modulated along the orbit, with a maximum close to periastron. We show that
the accretion itself cannot account for the periodical variation in
brightness. Instead, we explain the observed variation in the light curve with
absorption by the accretion tail. Our results are in general agreement with
previously derived stellar masses, orbital parameters, and the system
orientation, but imply that the B[e] supergiant is more luminous. We find an
effect related to the orbital motion of the two stars, in which the accretion
tail is cut by the primary and we term it the Lizard Autotomy Effect. As part
of the effect, the primary is self accreting wind that it ejected earlier. The
Lizard Autotomy Effect creates an outwardly expanding spiral shell made up of
broken segments. We suggest that such a tail exists in other B[e] supergiant
systems and can be the source of the circumstellar material observed in such
systems. The accretion also forms a disc around the secondary near periastron
that later vanishes. We suggest that the formation of such a disc can launch
jets that account for the bipolar structure observed around some B[e]
supergiants.

###### keywords:

stars: massive — stars: mass-loss — stars: winds, outflows — (stars:)
binaries: general — stars: emission-line, Be — accretion, accretion discs

††pubyear: 2023††pagerange: Accretion in the Binary System GG Carinae and
Implications for B[e] Supergiants–References

##  1 INTRODUCTION

Massive star evolution is an intriguing process that includes stages that are
not yet fully understood (e.g., Smartt, 2009; Owocki, 2011; Davidson &
Humphreys, 2012; Vink, 2015; Georgy et al., 2017; Owocki et al., 2019; Farrell
et al., 2022; Eldridge & Stanway, 2022). Amongst other massive stars at
evolved stages, the B[e] supergiant stars (B[e]SGs or sgB[e]) are rare stars,
and only a few tens of them are known. B[e]SGs exhibits optical spectra with
strong Balmer emission lines, low excitation permitted emission lines of
predominantly low ionization metals, and forbidden emission lines of [Fe II],
[O I] and other species, as well as strong infrared excess in their spectrum,
strong continua and few if any apparent absorption lines (e.g. Allen & Swings,
1976; Zickgraf et al., 1986, 1996; Lamers et al., 1998; Miroshnichenko, 2006;
Zickgraf, 2006; Aret et al., 2016; Humphreys et al., 2017; Oudmaijer &
Miroshnichenko, 2017; Kraus, 2019).

As such, and having approximately the same effective temperature but lower
luminosities ranging between $L=10^{4}$–$10^{6}~{}\rm{L_{\sun}}$, B[e]SGs may
be related as the fainter, lower mass siblings of Luminous Blue Variables
(LBVs) (e.g., Zickgraf et al., 1996; Clark et al., 2013; Kraus et al., 2014).
As the luminosity ranges overlap it may be difficult to distinguish between
them, with the main way being the lack of [O I] emission in LBVs (Humphreys et
al., 2017). Lower mass stars sometimes show similar properties to the B[e]
phenomenon, but not being SGs they do not belong to the group, and can be
classified as FS CMa stars (e.g., Miroshnichenko, 2006; Korčáková et al.,
2022; Miroshnichenko et al., 2023) and might be close binaries that are
undergoing or have undergone in the past a phase of mass transfer between the
components (see also Nodyarov et al. 2022 for a recent example of such a
system). Some yellow hypergiants might also appear as B[e]SGs when passing the
yellow void in the HR diagram (e.g., Davies et al., 2007; Aret et al., 2017a,
b).

B[e]SGs are surrounded with a circumstellar molecular and dusty disc or ring,
but the mechanism of how this circumstellar material arranged into its shape
and structure is not clear (Conti, 1997; Lamers et al., 1998; Kraus et al.,
2013; Kraus, 2019). The most stable and most abundant molecule emitting in the
near-IR in the discs around B[e]SGs is CO. The spectrum also shows emission
from the isotopic molecule 13CO, indicating that the circumstellar material
has been released from the evolved star (Kraus, 2009; Kraus et al., 2022). The
rarity of B[e]SGs may indicate that they are a fast transitional phase in
massive stars evolution of a certain mass range, characterized by strong mass
loss. A further strengthening to this conjecture comes from the fact that many
B[e]SGs show not only small-scale molecular and dusty rings, but also large-
scale nebulae and ejecta, suggesting previous mass-loss events (e.g., Kraus et
al., 2021).

B[e] stars have been found to also have large-scale (up to several pc) gas
with various shapes: spherical and ring nebulae, spiral-arms or partial spiral
arms, bipolar and unipolar-lobe structures (e.g., Marston & McCollum, 2008).
Radial velocities (RV) have been measured in the circumbinary material around
B[e] stars, suggesting an outflowing velocity component (Kraus et al., 2010).
Observations of the spatial distribution of the circumstellar gas revealed
that it is more likely accumulated in multiple rings or arcs, and possibly in
spiral arm-like structures (Torres et al., 2018; Kraus et al., 2016;
Maravelias et al., 2018; Liimets et al., 2022).

Some B[e]SGs were suggested to be in binary systems, and in this case the
circumstellar disc is in fact a circumbinary disc. According to Georgy et al.
(2017) it is difficult to explain the asymmetries in the circumstellar medium
of B[e]SGs with a single star, since their stellar evolution models predict
small rotation rates for blue supergiants at the observed position of B[e]SGs
in the H-R diagram. Thus, the similarity between the B[e]SGs, LBVs and FS Cma
might support a binary origin for the B[e]SGs as well.

GG Carinae (HD 94878) is one of these confirmed binary B[e]SG systems, that is
relatively well observed (Swings, 1974; Hernandez et al., 1981; Gosset et al.,
1985; Brandi et al., 1987; Machado et al., 2004; Pereyra et al., 2009;
Marchiano et al., 2012; Kraus et al., 2013; Porter et al., 2021, 2022). GG Car
is considered to be the archetypal B[e]SG system due to the high B[e]SG mass
$M_{1}=24\pm 4~{}\rm{M_{\sun}}$, luminosity
$\log(L_{1}/~{}\rm{L_{\sun}})=5.26\pm 0.19$ (Porter et al., 2021), and complex
circumstellar environment. However, its peculiar behavior on both short and
orbital time scales suggests that it is an atypical system. Most previous
evolutionary studies of GG Car assumed that the two binary components have
evolved as pseudo-single stars (Kraus, 2009, 2019; Kraus et al., 2013; Oksala
et al., 2013). Porter et al. (2021) brought a different view, modeling the
system as a binary with interaction between its components.

According to radiation transfer calculations with the cmfgen code computed by
Porter et al. (2021), the observed spectral behavior of GG Car implies that
the emission lines are generated by the primary’s wind, and the varying
amplitudes of the different line species are attributed to their formation at
different radii, with less energetic lines forming at larger radii. They also
determined the binary orbit by modeling the atmosphere of the primary and the
radial velocity variations of the wind lines. They concluded that the period,
which was earlier unclear, is $P=31.01\pm 0.01~{}\rm{days}$ and the
eccentricity is $e=0.50\pm 0.03$, much larger than earlier estimates
(Marchiano et al., 2012; Kraus et al., 2013). From their fit to the orbital
solution Porter et al. (2021) also found that the argument of periastron is
$\omega=339.87^{\circ}\pm 3.1^{\circ}$. Grant & Blundell (2022) used the
marginalized GP algorithm and found similar orbital eccentricity and
$\omega=334.36^{\circ}\pm 6.1^{\circ}$. The inclination of the system is
$i\approx 60^{\circ}$ (Kraus et al., 2013), but it is not very well
constrained, perhaps to within $20^{\circ}$ (Porter et al., 2021).

To explain the variation in luminosity along the orbit, Porter et al. (2021)
suggested that the system undergoes enhanced mass transfer and accretion from
the primary to the secondary at periastron. While the primary component of GG
Car is well identified, the class of the secondary component is still unknown,
though it is much less luminous than the primary. Porter et al. (2021)
obtained a preferred set of parameters for both stars in the system as well as
the binary orbit (values cited below).

The primary wind parameters obtained by Porter et al. (2021) are unusual as
the mass loss is relatively low and the wind is slow. Massive stars with
similar mass loss rate usually have higher wind velocities (e.g. Puls et al.,
2008; Owocki, 2015; Vink & Sander, 2021). We therefore test another set of
parameters based on the Geneva stellar evolution code (Ekström et al., 2012).
The simulations results that derived from using the two sets are then compared
to observations.

The GG Car system has a very weak x-ray emission, indicating that it is not a
colliding-wind binary such as $\eta$ Car (Pittard & Corcoran, 2002; Akashi et
al., 2006; Kashi & Soker, 2009a; Kashi et al., 2021) or a few other LBVs
(e.g., Nazé et al., 2012). Namely, the secondary is in an evolutionary state
in which its wind is very weak compared to the primary and has a low wind
momentum ratio. This indicates that the accretion is expected to be Bondi-
Hoyle-Lyttleton (BHL), or a fraction of this rate (e.g., Kashi et al., 2022).

The spectroscopic study of (Porter et al., 2022) found that GG Car has two
circumbinary rings. The rings are orbiting with projected velocities of 84.6
and 27.3 $~{}\rm{km~{}s^{-1}}$ at radii of $2.8^{+0.9}_{-1.1}~{}\rm{au}$ and
$27^{+9}_{-10}~{}\rm{au}$ for the inner and outer rings, respectively. They
further found variations in the rings consistent with pumping of the
eccentricity of a radially thin circumbinary ring by the inner binary. The
ring or rings are far from the binary and its immediate vicinity which are the
focus of our paper, but we will discuss how they might have formed.

In this paper we a run numerical simulations for the GG Car binary system
obtaining the geometry of the wind as a result of binary interaction as well
as accretion onto the secondary. In section 2 we describe the numerical
simulation. The results from our two sets of stellar and wind parameters are
presented in section 3. We discuss our results in section 4. Our summary and
conclusions are given in section 5.

##  2 THE NUMERICAL SIMULATIONS

Figure 1:  Left: Density maps with velocity vectors showing slices in the
orbital plane ($z=0$), at different times for set P. The times start after two
orbits in which we let the system relax and shows the third orbit from
periastron ($t=62~{}\rm{days}$), through apastron (middle right,
$t=77~{}\rm{days}$) till one day before the next periastron (lower right). The
periastron passage takes place when the primary is at
$(-0.305,0,0)~{}\rm{au}$, and moving counterclockwise when looking from
‘above’ the grid (positive $z$). The secondary is at the center, with its
radius indicated by a small black circle. The primary, with its radius
indicated by a large black circle, orbits the secondary counterclockwise, with
periastron occurring when the primary is left of the secondary. The accretion
tail is obtained and funnels the primary wind to be accreted onto the
secondary. Close to periastron passage ($t=62~{}\rm{days}$) the tail is being
cut by the passing of the secondary close to the primary. Right: Focusing on
the secondary at the same times as on the panels on the left hand side. Note
that the color scheme is the same for both panels but the velocity scale is
different, and more velocity vectors are plotted for the focus panels. An
accretion disc-like structure around the secondary is obtained shortly after
periastron passage. At later times the accretion is inhomogeneous BHL.

Figure 2:  Same as Figure 1, but for set G. The terminal velocity of the
primary for this set $v_{\rm 1,inf}=1400~{}\rm{km~{}s^{-1}}$ is much larger
than the orbital velocity $v_{\rm orb,p}<<v_{\rm 1,inf}$, Therefore a curved
tail that changes its direction along the orbits obtained. Differently from
model P, the primary does not cut the tail as in set P and the Lizard Autotomy
Effect does not take place.

We run the numerical simulation using the hydrodynamic code flash, described
in Fryxell et al. (2000). Our 3D Cartesian grid extends over
$x\in(-4,4)~{}\rm{au}$, $y\in(-4,4)~{}\rm{au}$ and $z\in(-2,2)~{}\rm{au}$. The
grid is centered around the secondary star, and the orbital plane is on the
$x$–$y$ plane at $z=0$. To solve the hydrodynamic equations we use the flash
version of the split piece-wise parabolic method (PPM) solver (Colella &
Woodward, 1984) which is an extension of the Godunov method, that solves a
Riemann problem and allows resolving shocks. A multigroup diffusion
approximation is used for solving the radiative transfer equation and updating
the energy. We also include radiative cooling from Sutherland & Dopita (1993),
which is necessary in modeling colliding winds governed by instabilities. The
temperature is limited from below to $1\,000~{}\rm{K}$.

We set up the GG Car binary system using two sets of stellar parameters. The
first set, designated ‘P’, is based on the parameters from in Porter et al.
(2021), noting that some of them were derived in earlier papers. The binary
semi-major axis is $a=0.61~{}\rm{au}$, the orbital period is $P=31.01$ days
and the eccentricity is $e=0.5$. The stellar masses of the primary B[e]SG star
and its companion are ${M_{1}=24~{}\rm{M_{\sun}}}$ and
${M_{2}=7.2~{}\rm{M_{\sun}}}$, respectively. We note that for the masses the
range in the literature was quite large. The effective temperature of the
primary is $T_{\rm eff,1}=23\,000~{}\rm{K}$, its radius is taken to be
$R_{1}=27~{}\rm{R_{\sun}}$ and the luminosity is $L_{1}\simeq 1.8\times
10^{5}~{}\rm{L_{\sun}}$. The secondary is hotter and smaller with $T_{\rm
eff,2}=16\,000~{}\rm{K}$, $R_{2}=7.2~{}\rm{R_{\sun}}$, and $L_{2}\simeq
3\,000~{}\rm{L_{\sun}}$. The composition used for both stars is appropriate
for an evolved star with CNO-processed material and the He abundance is taken
to be 50 percent, and the metalicity is $Z=0.02$. This composition is also
used for calculating the opacity, for which an opal model is used (Iglesias &
Rogers, 1996) and low temperature opacities are also included. The interiors
of the stars are excluded from the simulation. The gravitational field is
modeled by two point sources at the centers of both stars. The winds are
ejected from the grid cells in a narrow spherical shell around each star. Both
stars accelerate their winds according to the $\beta$ velocity law
$v(r)=v_{\infty}\left(1-R_{\ast}/r\right)^{\beta}$, where $v_{\infty}$ is the
terminal velocity, $R_{\ast}$ is the stellar radius, $r$ is the radial
distance measured from the center of the star and $\beta$ is the wind
acceleration parameter. The primary mass loss rate is $\dot{M}_{1}=2.2\times
10^{-6}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$ and its wind velocity has a terminal
value of $v_{1,\infty}=265~{}\rm{km~{}s^{-1}}$ with radiative acceleration
corresponding to $\beta=1$. The secondary mass loss rate is
$\dot{M}_{2}=10^{-8}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$, and its wind velocity
has a terminal value $v_{2,\infty}=1\,600~{}\rm{km~{}s^{-1}}$ with an
acceleration parameter $\beta=1$. The grid has 6 levels of refinement that are
kept constant with time. The largest cell has a side of $\simeq
18~{}\rm{R_{\sun}}$ and the smallest has a side of $\simeq
0.6~{}\rm{R_{\sun}}$. The periastron passage takes place when the primary is
at $(-0.305,0,0)~{}\rm{au}$, and moving counterclockwise when looking from
‘above’ the grid (positive $z$).

The second set of parameters we use, designated ‘G’, is based on a stellar
track from the Geneva code (Ekström et al., 2012). Using an evolutionary track
of a rotating star with zero age main sequence mass $M_{1,\rm
ZAMS}=25~{}\rm{M_{\sun}}$ for the primary, we set our simulation at
$M_{1}\simeq 23.6~{}\rm{M_{\sun}}$, when the effective temperature is $T_{\rm
1,eff}=23\,800~{}\rm{K}$, closest to the value in model P. The mass loss rate
is $\dot{M}_{1}\simeq 4\times 10^{-6}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$, and
the terminal velocity we derive is $v_{1,\infty}=2\simeq
1400~{}\rm{km~{}s^{-1}}$. The main difference between the two sets of
parameters is therefore that $v_{1,\infty}$ is much larger for set G. For the
secondary we use $M_{2,\rm ZAMS}=7.2~{}\rm{M_{\sun}}$ which gives almost the
same mass at the time we take for the primary, as it is still a main-sequence
star, $M_{2}\simeq 7.2~{}\rm{M_{\sun}}$. The terminal wind velocity is
$v_{2,\infty}\simeq 900~{}\rm{km~{}s^{-1}}$ and the mass loss rate is
$\dot{M}_{2}\simeq 6.24\times 10^{-12}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$. Table
1 gives the full list of parameters in both models. The orbital parameters of
set G are similar to those of set P (except for a minor difference of $<0.7$%
in the semi-major axis). As the secondary is smaller, for set G the grid has 8
levels of refinement, and the smallest has a side of $\simeq
0.14~{}\rm{R_{\sun}}$.

Table 1:  List of stellar, wind and and orbital parameters for the two sets of parameters used in our simulations. Set P is based on parameters from Porter et al. (2021). Set G is based on a stellar track from the Geneva code.  Parameter | Set P | Set G
---|---|---
$M_{1}(~{}\rm{M_{\sun}})$ | $24$ | $23.6$
$R_{1}(~{}\rm{R_{\sun}})$ | $27$ | $31$
$T_{\rm eff,1}(~{}\rm{K})$ | $23\,000$ | $23\,800$
$L_{1}(~{}\rm{L_{\sun}})$ | $1.8\times 10^{5}$ | $2.4\times 10^{5}$
$\dot{M}_{1}(~{}\rm{M_{\odot}}~{}\rm{yr^{-1}})$ | $2.2\times 10^{-6}$ | $4\times 10^{-6}$
$v_{1,\infty}(~{}\rm{km~{}s^{-1}})$ | $265$ | $1\,400$
$\beta_{1}$ | $1$ | $1$
$M_{2}(~{}\rm{M_{\sun}})$ | $7.2$ | $7.2$
$R_{2}(~{}\rm{M_{\sun}})$ | $7.2$ | $3.4$
$T_{\rm eff,2}(~{}\rm{K})$ | $16\,000$ | $20\,800$
$\dot{M}_{2}(~{}\rm{M_{\odot}}~{}\rm{yr^{-1}})$ | $10^{-8}$ | $6.24\times 10^{-12}$
$v_{2,\infty}(~{}\rm{km~{}s^{-1}})$ | $1\,600$ | $900$
$\beta_{2}$ | $1$ | $1$
$P$ (days) | $31.01$ | $31.01$
$e$ | $0.5$ | $0.5$
$a(~{}\rm{au})$ | $0.61$ | $0.61$

For both sets of parameter, the wind parameters of the secondary have a very
small impact on the flow structure of the simulation as its wind momentum is
negligible compared to that of the primary. So the grid is dominated by almost
exclusively the primary wind, that is shaped by the orbital motion and the
gravity of the secondary. We ran one simulation for each set of parameters
with one set of parameters, but post-process it for four different
orientations.

Our simulation code does not assume that accretion takes place. We design the
code to allow treatment of accretion if material flows towards the secondary
and not stopped by radiation braking. The accretion treatment includes
measuring the accreted material, its arrival direction and other properties.
We quantify the accretion rate according to the method described in (Kashi,
2019). The secondary is resolved as a sphere (not a sink) and accretion is
calculated in each cell around the secondary.

##  3 RESULTS

###  3.1 Set P

Figure 3:  A 3D view of the simulation with parameters set P. The accretion
tail being cut by the primary during periastron passage, shown from two
viewing angles. The secondary does not appear in the 3D view, and was added
for perspective and to show the scale ($R_{2}=7.2~{}\rm{R_{\sun}}$). There is
a dense tail behind the secondary, part of which is an accretion tail towards
the secondary and part is outflowing away from it. The secondary takes a sharp
elliptical turn close to periastron and as a result the close primary and its
wind cut the tail and creates a gap in the tail – the Lizard Autotomy Effect.
The process repeats every periastron passage, creating an outflowing spiral
with gaps.

We start with set P. The simulation outputs are exported and post-processed
every 0.5 day. Figure 1 shows a slice of the simulation grid on the orbital
plane. The right six panels are zoom-in on the secondary at times
corresponding to the left six panels, and have the same density scale, but
more velocity vectors. The times given in Figure 1 in the upper left of each
panel are during the third orbit, as for two orbital periods we let the system
relax from the initial condition that did not include wind from previous
orbits but only wind that started to flow from the two stars. Therefore after
two orbits the wind structure is representative of the system unlike the
initialization of the simulation. The times show the third orbit from
periastron (upper left panel in each six panel group), through apastron
(middle right) till one day before the next periastron (lower right).

As expected, we find that the accretion takes a geometric form that resembles
inhomogeneous BHL accretion (Livio et al., 1986). The primary wind flows from
both sides of the secondary and focused by its gravity into a curved tail.
Part of the accretion tail flows towards the secondary and gets accreted, and
part of it flows outwards and creates a spiral structure. Accretion also takes
place from other directions but at a lower rate. After periastron passage a
clear dense accretion disc is obtained around the secondary (see the middle
left panel at $t=65.0~{}\rm{days}$). The disc lasts for a few more days and
dissipates as the system approaches apastron (middle right panel at
$t=77.0~{}\rm{days}$). An accretion tail remains even after the dissipation of
the disc, and lasts until the next periastron passage. The flow creates a
spiral, with higher density close to the secondary.

Figure 3 shows a 3D contour plot of the high density parts of the simulation
close to periastron passage. These dense parts include the strong wind of the
primary and the tail of the secondary. There is a dense tail behind the
secondary, part of which is an accretion tail towards the secondary and part
is outflowing away from it. We can see the secondary (added artificially on
the right side that views the system from above the orbital plane) close to
periastron passage. As the secondary takes a sharp elliptical turn when coming
to periastron and passes close to the primary and its wind, the tail is being
cut by the primary and a gap in the tail is created. As this process is
reminiscent of the self cutting of a lizard’s tail we will refer to it as the
Lizard Autotomy Effect. The process repeats every periastron passage, creating
an outflowing spiral with gaps, that also has a rotating velocity component.

Another effect that we obtain as part of the Lizard Autotomy Effect is self
accretion: the primary passes through the tail and accretes its own wind. We
did not measure quantitatively the accreted mass onto the primary, as our code
is not designed to do so (it only measures the mass accretion onto the
secondary). The self accreted material comes mostly from one direction – the
side of the primary facing the tail. As the density of the tail is much larger
than the density of the direct primary wind (Figure 3) it is unlikely that the
latter can prevent the accretion.

Figure 4 shows the accretion rate onto the secondary as a function of time.
The accretion we found is modulated by the orbital period. Accretion occurs
throughout the binary period but much stronger close to periastron. When the
system is close to apastron there is still accretion taking place from
residual material that the primary left around the secondary. There is also a
temporary formation of an accretion disc around the secondary that dissipates
shortly after periastron passage. The residual accertion is from filaments of
gas compressed by the orbital motion that remain gravitationaly bound to the
secondary.

Figure 4:  The mass accretion rate obtained in our simulation or set P. The
first two orbits ($t<62~{}\rm{days}$) are the the initialization of the
simulation, where the structure of the wind does not yet represent the system.
We ran the simulation for seven orbitals period (the first three are shown in
the figure) and from the third orbit we obtain repetitive behaviour of the
accretion rate from orbit to orbit. We can see the accretion rate is periodic
with the orbit, with a peak just before periastron. A sharp decline following
the peak indicates the effect we term the Lizard Autotomy Effect, in which the
primary cuts the accretion tail, and by that reduces the mass accretion rate.
Figure 5: Similar to Figure 4, but for the G parameters set. The accretion
rate onto the secondary is $\propto r^{-2}$. The Lizard Autotomy Effect is not
evident. The mass accretion rate is much smaller that for set P.

Porter et al. (2021) suggested that the accretion of the primary wind by the
secondary is the cause of the photometric variations at the orbital period in
GG Car. The photometric variability according to their observations is $\Delta
m_{V}\simeq 0.2\rm{mag}$, which corresponds to a decrease of about 20% in
luminosity. This conclusion was reached after excluding other options like
extinction, scattering, and reflection by circumbinary material. While
accretion onto the companion can generally account for increase in brightening
of a star or dimming of this light, we argue that in the case of GG Car
accretion by itself cannot cause the observed brightening.

Using the mass loss rate of the primary $\dot{M}_{1}=2.2\times
10^{-6}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$ from Porter et al. (2021), the
resulted highest mass accretion rate onto the secondary obtained by our
simulations is $\dot{M}_{\rm 2,acc}\simeq 8\times
10^{-7}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$. The resulted accretion luminosity is

$\begin{split}L_{\rm 2,acc}&\simeq
25\left(\frac{M_{2}}{7.2~{}\rm{M_{\sun}}}\right)\left(\frac{\dot{M}_{\rm
2,acc}}{8\times 10^{-7}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}}\right)\\\
&\times\left(\frac{R_{2}}{7.2~{}\rm{R_{\sun}}}\right)^{-1}~{}\rm{L_{\sun}},\end{split}$
(1)

which is a negligible fraction of the primary luminosity $L_{\rm
2,acc}/L_{1}\simeq 1.4\times 10^{-4}$. Moreover, even if all the mass lost
from the primary gets accreted onto the secondary it will still not be enough
to explain an increase of $\simeq 20$% in luminosity. The other option,
dimming of the secondary by the accreted material, is also not enough to
account for the observed increase in luminosity, as the secondary luminosity
is at most $4800~{}\rm{L_{\sun}}$, only $<3\%$ of the primary luminosity, and
even if obscured completely will not suffice.

We therefore suggest another explanation – attenuation of the primary flux by
the surrounding material. We suggest that the luminosity does not periodically
increase as a result of the accretion onto the companion, but rather
periodically decreases as a result of absorption by the surrounding material,
mainly in the tail. The origin of the surrounding material is the wind of the
primary itself. The gravity of the companion and the eccentric orbital motion
creates an accretion structure with disc and tail, as shown in Figure 3. We
take the emission to originate at the surface of the primary star, and
calculate the optical depth along the line of sight till reaching a boundary
of the simulation box

$\tau(t)=\int\kappa(t)\rho(t)\,dl.$ (2)

We then calculate the attenuation of the observed magnitude by

$M_{V}(t)=M_{V,0}+A_{V}(t)=M_{V,0}+1.086\tau(t),$ (3)

Where $M_{V,0}$ is the non-obscured magnitude of the star and $M_{V}$ is the
observed magnitude through the absorbing material.

Figure 6:  The simulation results for the light curve of GG Car for parameters
set G. The observed magnitude is taken from OMC (blue circles) and ASAS (blue
crosses), and folded to the time of the orbit according to the ephemeris in
Porter et al. (2021). The observation results are binned into multiple bins,
each includes 1/15 of the orbit, symmetric in both sides. For each bin the
average and the standard deviation are calculated. The error bars indicate
$1\sigma$. The light curve follows the binary orbit. The tail and the spiral
structure of the wind attenuate the light of the primary. The secondary is a
negligible source. The sharp rise close to periastron passage (phases 0, 1, 2)
is the result of the primary cutting the tail and being exposed to the
observer. The resulted light curve from our simulation is given for different
inclination angles and $\omega$. We find the orientation $i=60^{\circ}$ and
$\omega=340^{\circ}$ gives the best results.

The results of our simulation for set P on top of the observed magnitude is
shown in Figure 6. The observed visible magnitude is taken from two sources
and combined together, archival observations from the Optical Monitoring
Camera (OMC) on-board the INTEGRAL satellite (Domingo et al., 2010), marked
with blue circles and The All Sky Automated Survey (ASAS) catalog (Pojmanski,
2003), marked with blue crosses. The observation results are binned into
multiple bins, using a running window. Each bin includes 1/15 of the orbit,
symmetric in both sides. For each bin the average and the standard deviation
are calculated. The result is shown in blue squares with 1$\sigma$ error-bars.
The resulted light curve from our simulation is over plotted. We used four
viewing angles with different inclination angles and $\omega$, within the
range obtained by Porter et al. (2021) and Grant & Blundell (2022). We find
the $i=60^{\circ}$ and $\omega=340^{\circ}$ gives the best results out of the
four orientations we attempted, in agreement with (Porter et al., 2021).

The implication of our results is that the primary luminosity is in fact
$\simeq 25\%$ larger than previously derived, namely
$\log(L_{1}/~{}\rm{L_{\sun}})=5.36\pm 0.2$. This also indicates that the mass
of the primary of GG Car is larger, and by extension also that of the
secondary, as the derivation of the latter relies on the mass ratio between
the two stars Porter et al. (2021).

The tail and the spiral structure obtained in our simulations should also be
evident in other observations. To show that, we process our simulation for set
P to obtain the RV of the absorption component for the He I lines. We assume
that the lines originate on the surface of the primary. For each length $dl$
along the line of sight to the primary we calculate $d\tau$ and the projected
velocity of the gas along the path towards the observer
$\mathbf{v_{p}}=(\mathbf{v}\cdot\mathbf{dl})d\hat{l}$, where $d\hat{l}$ points
to the line of sight, namely $i=60^{\circ}$ and $\omega=340^{\circ}$. We thus
obtain $d\tau(v_{p})$ along the line of sight, and get the radial velocity of
the absorption. We repeat the process for every $0.5$ day along the orbital
period.

Porter et al. (2021) showed three trailed spectra of He I lines along the
orbital period. The absorption RV showed periodic variation, having larger
values close to periastron. In Figure 7 we plot our simulation result over
figure 4 of Porter et al. (2021) that shows the trailed spectra of He I 6678
taken by the Global Jet Watch (GJW) telescopes. As we do not use advanced
radiation transfer procedure a perfect match to a specific line is not
expected. However, our results show that the RV obtained from the simulation
(brown line) generally follows the same behaviour, having minimal values close
to apastron and maximal values close to periastron passage.

Figure 7:  The result of our simulation for the RV of the absorped projected
velocity, assuming line of sight in the direction $i=60^{\circ}$ and
$\omega=340^{\circ}$ (brown line). The results are plotted on top of the
radial velocity profile of the He I 6678 line from Porter et al. (2021), taken
by GJW (the color background with bluer colors represent absorption and
yellower colors represent emission; precise colormaps were not provided). The
simulation results follow the observed absorption component, showing minimum
values close to apastron (phase 0.5, 1.5) and maximal values close to
periastron passage (phase 0, 1, 2).

###  3.2 Set G

We repeat the analysis for parameters set G. Figure 2 shows the density plots
for set G. For a large terminal velocity $v_{\rm
1,inf}=1400~{}\rm{km~{}s^{-1}}$ (set G) the orbital velocity at periastron is
much smaller $v_{\rm orb,p}<<v_{\rm 1,inf}$. Therefore a curved tail is
obtained but the primary does not cut as in set P and the Lizard Autotomy
Effect is non existent. The result is a flow with a narrow curved tail behind
the secondary, that changes relatively slowly along the orbit. To have a
similar Lizard Autotomy Effect as in set P, a much larger eccentricity is
required.

Therefore for set G the mass accretion rate onto the secondary changes almost
as $r^{-2}$ due to the orbital separation and no other effects are involved.
The obtained accretion rate is shown in Figure 5. The peak accretion rate for
set G is $\dot{M}_{\rm 2,acc}\simeq 3.4\times
10^{-8}~{}\rm{M_{\odot}}~{}\rm{yr^{-1}}$, much smaller than for set P. The
reason is the absence of the accretion tail in set G, that leads to no
efficient focusing mechanism of material into the secondary.

In Figure 8 we calculate the obscuration from different lines of sight. We
obtain that the magnitude of variation is much smaller than the observed one,
despite the mass loss rate of the primary being almost twice larger than in
set P. No result from the lines of sights that we checked (in any inclination)
was found to reproduce the observed periodical light curve. Even for
$i=90^{\circ}$, where a maximal eclipse is obtained, the decline is not enough
to account for the observed light curve. Moreover, even if we were to change
the line of sight such that the eclipse is at phase 0.5 it would not be
enough, as the observed decline is larger.

Figure 8:  The simulation results for the light curve of GG Car based on
parameters set G. The observations are as in Figure 6. The resulted light
curve from our simulation is given for different inclination angles and
$\omega$. The obscuration is periodic and follow the orbit. For $i=90^{\circ}$
an eclipse is seen at phase $\simeq 0.87$. We find that for none of the lines
of sight the obscuration is sufficient to account for the decline in the light
curve.

##  4 DISCUSSION

Between the two sets of stellar and wind parameters that we tested, we found
that the set of parameters we adopted from Porter et al. (2021) is a much
better fit with the observations. For this set we found that the primary cuts
the accretion tail close to periastron passage – a process we termed the
Lizard Autotomy Effect. A representation of this effect is shown in Figure 3.
The Lizard Autotomy Effect that we obtained in this paper causes the primary
to self accrete mass from wind that it had blown earlier and formed the tail,
and creates an expanding spiral shell with broken segments that expands
outwards.

We suggest that such tails exists in other B[e]SGs systems and can be the
source of the circumstellar material observed in such systems. The accretion
onto the secondary also forms a disc around the secondary (not a circumbinary
disc) near periastron that later disperses. Some of the gas later remains, and
form dense filament (as seen in Figure 1). We suggest that the formation of
such a disc can be accompanied by the launching of two opposite jets
perpendicular to the plane of the disc. These two jets can push ambient gas
and form a bipolar structure.

Such bipolar structures were observed around some B[e]SGs. For example,
Liimets et al. (2022) showed an H$\alpha$ image of the nebula around MWC 314,
a well known B[e]SG, that clearly shows its bipolar structure. Marston &
McCollum (2008) observed nebular features of the B[e] star MWC 819 seen in
2001, later to be confirmed as bipolar in a deeper image taken in 2020 by
Liimets et al. (2022). The B[e] star LHA 120-S 35 exhibits noticeable changes
in its spectral patterns in both the optical and infrared regions and its
P-Cygni line profiles of H I, Fe II, and O I indicate the presence of a robust
bipolar clumped wind (Torres et al., 2018). The B[e]SG MWC 137 is surrounded
by an optical nebula that was modeled as consisting of two double cones,
though might be related to stellar pulsations (Kraus et al., 2017, 2021).
Sánchez Contreras et al. (2019) suggested an accretion disc around a compact
object and jet-launching as a model for forming the red square nebula around
the B[e] star MWC 922.

It is well established that jets can shape bipolar nebulae. The process
applies to various astrophysical objects (Frank, 1999) including planetary
nebulae (e.g., Soker, 1990, 1998, 2002; Soker & Livio, 1994; Frank et al.,
1996; Sahai & Trauger, 1998; Akashi & Soker, 2008; Soker & Kashi, 2012; Balick
et al., 2019), post-AGB nebulae that involve a sudden anisotropic ejection of
stellar gas by a very late AGB (or very early post-AGB) star during a common-
envelope phase (e.g., Lee et al., 2009; Sahai et al., 2017; Soker, 2020; Zou
et al., 2020; Bollen et al., 2022; López-Cámara et al., 2022), young stellar
objects (e.g., Shu et al., 1994; Ray & Ferreira, 2021) Luminous Blue Variables
(e.g., Nota et al., 1995; Soker, 2001; Kashi & Soker, 2010), core collapse
supernovae (e.g., Gilkis et al., 2019; Akashi & Soker, 2021; Soker, 2022;
Akashi et al., 2023), and intermediate-luminosity optical transients (e.g.,
Soker & Kashi, 2016). For planetary nebulae there are also indications that
impulsive ejection of gas creates a bipolar structure (e.g., Akashi & Soker,
2008), a similar idea to what we propose here for GG Car and similar B[e]SGs.

The evidence for bipolar structure around B[e]SGs implies that very likely
jets ejected from these stars were involved in their formation. This, in turn,
implies there should exist an accretion disc to form these jet that exists
either continuously or intermittently. Our simulation obtained such an
accretion disc around the secondary that exists for every orbit close to
periastron passage, and therefore support the observed bipolar structure of
nebulae around B[e]SGs.

As earlier noted, some B[e]SGs show large-scale nebulae and ejecta, suggesting
previous mass-loss events (e.g., Kraus et al., 2021). These might be the
result of the expanding spiral arms that dilute and merge with those of
previous cycles. The spectroscopic study of Porter et al. (2022) found that GG
Car has two circumbinary rings. These rings are located further out from the
simulation grid. However, the outflow we find, in the shape of a spiral moving
outward, could have formed such rings further out. This is likely if the
spirals decelerate and accumulate in the radial direction, to form ring
structures.

We suggest that the circumstellar disc around B[e] stars, or at least some
B[e] stars, is in fact a cropped spiral created by the interaction of the
stellar wind with the companion star. The results here demonstrate that
accretion is an important process in massive binaries, and should be
considered in every close binary system that has uneven stellar wind momentum,
and especially in ones with eccentric orbit.

##  5 SUMMARY and CONCLUSIONS

We used hydrodynamic simulations to show that the B[e]SG binary system GG Car
experiences accretion onto the secondary close to periastron passage. By that
we confirm the earlier suggestion for this process to occur in this system by
Porter et al. (2021). We calculated the accretion rate onto the secondary star
and found peaks in the accretion rate (Figure 4) as suggested by Porter et al.
(2021). An accretion disc-like structure around the secondary is obtained
shortly after periastron passage, and dissipates as the system approaches
apastron (Figure 1). We also found that an asymmetric spiral structure is
formed around the system.

We showed that the variation in the V light curve is not as a result of an
increase to the luminosity from accretion, as the accretion power is not
sufficient to account for it. Instead we propose that the primary is more
luminous and it is dimmed periodically as the mass that it lost is compressed
into the accretion tail and the spiral structure which adsorbs its light
(Figure 6). The implication is that the primary of GG Car is more luminous by
$\approx 25\%$ than previously derived, and consequently that the stellar
masses are larger.

Accretion close to periastron passage is known to occur in other stars. One
famous example is $\eta$ Car, an LBV star that shares similar properties with
GG Car. Close to periastron in the eccentric ($e=0.9$) and long
($P=5.54~{}\rm{yrs}$) orbit the star accretes mass from the equatorial region
(Kashi & Soker, 2008; Kashi, 2017), and a disc is formed and stays for a few
weeks, and accounts for multiple observations across the spectra (Kashi &
Soker, 2007a, 2009b), including the variation in He lines (Kashi & Soker,
2007b), as in GG Car.

Porter et al. (2021) suggested that the fact that the blueshifted absorption
component in GG Car becomes deeper around periastron passage likely indicates
that more mass is being lost from the primary’s wind at that time. They
concluded that the primary mass-loss rate and the rate of mass being captured
by the secondary strongly depends on the phase of the binary due to the
eccentric orbit, and therefore we can expect that there will be changes in the
brightness due to varying rates of mass transfer and accretion. In our
simulation we kept the mass loss rate of the primary fixed, and did not
enhance it near periastron. We agree with Porter et al. (2021) that the mass
loss enhancement is reasonable. However, the mass accretion rate is modulated
according to the orbit regardless of the primary mass loss rate, and this is
indeed obtained for the two sets of parameters we used. Enhancing the primary
mass loss rate close to periaston can let alone enhance the accretion, but
accretion exists even without the enhancement.

In the present simulation we did not modify the stellar atmosphere of either
star as a result of accretion. Accretion at a very high rate can lead to a
reduction in the effective temperature of the star and in turn weaken its
wind, as suggested by Kashi & Soker (2009b) for accretion close to the
periastron passage of $\eta$ Carinae. This effect was later obtained in the
simulations of Kashi (2017).

In this paper we tried two representative sets of parameters. We did not make
an attempt to run the simulation with tuning the parameters, such as stellar
masses, radii, effective temperatures and eccentricity. We showed that one of
these parameter set (set P) reasonably reproduce photometric and spectroscopic
observations. The parameter set that reproduced the observations has a lower
primary star wind velocity then expected from stellar evolution models,
revealing an interesting problem that is beyond the scope of this work.

There are many spectroscopic observations for GG Car, and our simulation and
following analysis is limited and cannot reproduce all of them. However, we
showed that the RV absorption of He I 6678 matches the results of our
simulation (Figure 7). This line is emitted from the primary or its wind, very
close to the stellar surface (Porter et al., 2021), depends on the line of
sight and less sensitive to other parameters. Running more simulations on a
set of parameters can very likely improve the constraints on the stellar and
binary parameters. Likewise, processing the results with radiation transfer
code can give much better reproduction of spectral line observations than the
relatively simple procedures we performed in this work.

Our simulations show that binary interaction can have a significant role in
the formation of B[e]SGs and in shaping their circumstellar ans circumbinary
environment. Hydrodynamic simulations are an important tool for analysing
massive binary systems with complex wind structure, and in particular B[e]SGs,
as done in this case.

## Acknowledgements

We thank an anonymous referee for helpful comments. We acknowledge support
from the R&D Authority at Ariel University. We acknowledge the Ariel HPC
Center at Ariel University for providing computing resources that have
contributed to the research results reported within this paper.

## Data Availability

OMC V-band photometry observations are available at https://sdc.cab.inta-
csic.es/omc/secure/formbusqueda.jsp.
ASAS V-band photometry observations are available at
http://www.astrouw.edu.pl/cgi-asas/asascgigetdata?105559-6023.5,asas3.
This paper used Geneva code stellar evolution tracks from SYCLIST
https://www.unige.ch/sciences/astro/evolution/en/database/syclist/.
The data underlying this article will be shared on a reasonable request from
the corresponding author.

## References

  * Akashi et al. (2023) Akashi, M., Michaelis, A., & Soker, N. 2023, MNRAS, 518, 6123
  * Akashi & Soker (2008) Akashi, M., & Soker, N. 2008, MNRAS, 391, 1063
  * Akashi & Soker (2021) Akashi, M., & Soker, N. 2021, MNRAS, 501, 4053
  * Akashi et al. (2006) Akashi, M., Soker, N., & Behar, E. 2006, ApJ, 644, 451
  * Allen & Swings (1976) Allen, D. A., & Swings, J. P. 1976, A&A, 47, 293
  * Aret et al. (2016) Aret, A., Kraus, M., & Šlechta, M. 2016, MNRAS, 456, 1424
  * Aret et al. (2017a) Aret, A., Kraus, M., Kolka, I., & Maravelias, G. 2017, Stars: From Collapse to Collapse, 510, 162
  * Aret et al. (2017b) Aret, A., Kolka, I., Kraus, M., & Maravelias, G. 2017, The B[e] Phenomenon: Forty Years of Studies, 508, 239
  * Balick et al. (2019) Balick, B., Frank, A., & Liu, B. 2019, ApJ, 877, 30
  * Bollen et al. (2022) Bollen, D., Kamath, D., Van Winckel, H., et al. 2022, A&A, 666, A40
  * Brandi et al. (1987) Brandi, E., Gosset, E., & Swings, J.-P. 1987, A&A, 175, 151
  * Clark et al. (2013) Clark, J. S., Bartlett, E. S., Coe, M. J., et al. 2013, A&A, 560, A10
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174
  * Conti (1997) Conti, P. S. 1997, Luminous Blue Variables: Massive Stars in Transition, 120, 161
  * Davidson & Humphreys (2012) Davidson, K., & Humphreys, R. M. 2012, Eta Carinae and the Supernova Impostors, 384.
  * Davies et al. (2007) Davies, B., Oudmaijer, R. D., & Sahu, K. C. 2007, ApJ, 671, 2059
  * Domingo et al. (2010) Domingo, A., Gutiérrez-Sánchez, R., Rísquez, D., et al. 2010, Highlights of Spanish Astrophysics V, 14, 493
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
  * Eldridge & Stanway (2022) Eldridge, J. J., & Stanway, E. R. 2022, ARA&A, 60, 455
  * Farrell et al. (2022) Farrell, E., Groh, J. H., Meynet, G., & Eldridge, J. J. 2022, MNRAS, 512, 4116.
  * Frank (1999) Frank, A. 1999, New Astron. Rev., 43, 31
  * Frank et al. (1996) Frank, A., Balick, B., & Livio, M. 1996, ApJ, 471, L53
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
  * Georgy et al. (2017) Georgy, C., Saio, H., Ekström, S., & Meynet, G. 2017, The B[e] Phenomenon: Forty Years of Studies, 508, 99
  * Gilkis et al. (2019) Gilkis, A., Soker, N., & Kashi, A. 2019, MNRAS, 482, 4233
  * Glatzel & Mehren (1996) Glatzel, W., & Mehren, S. 1996, MNRAS, 282, 1470
  * Gosset et al. (1985) Gosset, E., Hutsemekers, D., Surdej, J., & Swings, J. P. 1985, A&A, 153, 71
  * Grant & Blundell (2022) Grant, D., & Blundell, K. 2022, MNRAS, 509, 367
  * Hernandez et al. (1981) Hernandez, C. A., Lopez, L., Sahade, J., & Thackeray, A. D. 1981, PASP, 93, 747
  * Humphreys et al. (2017) Humphreys, R. M., Gordon, M. S., Martin, J. C., Weis, K., & Hahn, D. 2017, ApJ, 836, 64
  * Iglesias & Rogers (1996) Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943.
  * Kashi (2017) Kashi, A. 2017, MNRAS, 464, 775
  * Kashi (2019) Kashi, A. 2019, MNRAS, 486, 926
  * Kashi et al. (2021) Kashi, A., Principe, D. A., Soker, N., & Kastner, J. H. 2021, ApJ, 914, 47
  * Kashi & Soker (2007a) Kashi, A., & Soker, N. 2007, MNRAS, 378, 1609
  * Kashi & Soker (2007b) Kashi, A., & Soker, N. 2007, New Astron., 12, 590
  * Kashi & Soker (2008) Kashi, A., & Soker, N. 2008, New Astron., 13, 569
  * Kashi & Soker (2009a) Kashi, A., & Soker, N. 2009, MNRAS, 397, 1426
  * Kashi & Soker (2009b) Kashi, A., & Soker, N. 2009, New Astron., 14, 11
  * Kashi & Soker (2010) Kashi, A., & Soker, N. 2010, ApJ, 723, 602
  * Kashi et al. (2022) Kashi, A., Michaelis, A., & Kaminetsky, Y. 2022, MNRAS, 516, 3193
  * Korčáková et al. (2022) Korčáková, D., Sestito, F., Manset, N., et al. 2022, A&A, 659, A35
  * Kraus (2009) Kraus, M. 2009, A&A, 494, 253
  * Kraus (2019) Kraus, M. 2019, Galaxies, 7, 83
  * Kraus et al. (2010) Kraus, M., Borges Fernandes, M., & de Araujo, F. X. 2010, A&A, 517, A30
  * Kraus et al. (2016) Kraus, M., Cidale, L. S., Arias, M. L., et al. 2016, A&A, 593, A112
  * Kraus et al. (2014) Kraus, M., Cidale, L. S., Arias, M. L., Oksala, M. E., & Borges Fernandes, M. 2014, ApJ, 780, L10
  * Kraus et al. (2022) Kraus, M., Cidale, L. S., Arias, M. L., et al. 2022, IAU Symposium, 366, 51
  * Kraus et al. (2017) Kraus, M., Liimets, T., Cappa, C. E., et al. 2017, AJ, 154, 186
  * Kraus et al. (2021) Kraus, M., Liimets, T., Moiseev, A., et al. 2021, AJ, 162, 150
  * Kraus et al. (2013) Kraus, M., Oksala, M. E., Nickeler, D. H., et al. 2013, A&A, 549, A28
  * Lamers et al. (1998) Lamers, H. J. G. L. M., Zickgraf, F.-J., de Winter, D., Houziaux, L., & Zorec, J. 1998, A&A, 340, 117
  * Lee et al. (2009) Lee, C.-F., Hsu, M.-C., & Sahai, R. 2009, ApJ, 696, 1630
  * Liimets et al. (2022) Liimets, T., Kraus, M., Moiseev, A., et al. 2022, Galaxies, 10, 41
  * Livio et al. (1986) Livio, M., Soker, N., de Kool, M., & Savonije, G. J. 1986, MNRAS, 222, 235
  * López-Cámara et al. (2022) López-Cámara, D., De Colle, F., Moreno Méndez, E., Shiber, S., & Iaconi, R. 2022, MNRAS, 513, 3634
  * Machado et al. (2004) Machado, M. A., de Araújo, F. X., de Faria Lopes, D., & Pereira, C. B. 2004, Revista Mexicana de Astronomia y Astrofisica Conference Series, 20, 239
  * Marchiano et al. (2012) Marchiano, P., Brandi, E., Muratore, M. F., et al. 2012, A&A, 540, A91
  * Maravelias et al. (2018) Maravelias, G., Kraus, M., Cidale, L. S., et al. 2018, MNRAS, 480, 320
  * Marston & McCollum (2008) Marston, A. P., & McCollum, B. 2008, A&A, 477, 193
  * Miroshnichenko (2006) Miroshnichenko, A. S. 2006, Stars with the B[e] Phenomenon, 355, 13
  * Miroshnichenko et al. (2023) Miroshnichenko, A.S.; Zharikov, S.V.; Manset, N.; Khokhlov, S.A.; Nodyarov, A.S.; Klochkova, V.G.; Danford, S.; Kuratova, A.K.; Mennickent, R.; Chojnowski, S.D.; Raj, A.; Bisht, D. Galaxies 2023, 11, 36.
  * Nazé et al. (2012) Nazé, Y., Rauw, G., & Hutsemékers, D. 2012, A&A, 538, A47
  * Nodyarov et al. (2022) Nodyarov, A. S., Miroshnichenko, A. S., Khokhlov, S. A., et al. 2022, ApJ, 936, 129
  * Nota et al. (1995) Nota, A., Livio, M., Clampin, M., & Schulte-Ladbeck, R. 1995, ApJ, 448, 788
  * Oksala et al. (2013) Oksala, M. E., Kraus, M., Cidale, L. S., Muratore, M. F., & Borges Fernandes, M. 2013, A&A, 558, A17
  * Oudmaijer & Miroshnichenko (2017) Oudmaijer, R. D., & Miroshnichenko, A. S. 2017, The B[e] Phenomenon: Forty Years of Studies, 508, 3
  * Owocki (2011) Owocki, S. 2011, Bulletin de la Societe Royale des Sciences de Liege, 80, 16
  * Owocki (2013) Owocki, S. 2013, Planets, Stars and Stellar Systems. Volume 4: Stellar Structure and Evolution, 4, 735 :Stellar Winds
  * Owocki (2015) Owocki, S. P. 2015, Very Massive Stars in the Local Universe, 412, 113
  * Owocki et al. (2019) Owocki, S. P., Hirai, R., Podsiadlowski, P., & Schneider, F. R. N. 2019, MNRAS, 485, 988
  * Pereyra et al. (2009) Pereyra, A., de Araújo, F. X., Magalhães, A. M., Borges Fernandes, M., & Domiciano de Souza, A. 2009, A&A, 508, 1337
  * Pittard & Corcoran (2002) Pittard, J. M., & Corcoran, M. F. 2002, A&A, 383, 636
  * Pojmanski (2003) Pojmanski, G. 2003, Acta Astron., 53, 341
  * Porter et al. (2021) Porter, A., Grant, D., Blundell, K., & Lee, S. 2021, MNRAS, 501, 5554
  * Porter et al. (2022) Porter, A., Blundell, K., & Lee, S. 2022, MNRAS, 509, 1720
  * Puls et al. (2008) Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209
  * Ray & Ferreira (2021) Ray, T. P., & Ferreira, J. 2021, New Astron. Rev., 93, 101615
  * Sahai & Trauger (1998) Sahai, R., & Trauger, J. T. 1998, AJ, 116, 1357
  * Sahai et al. (2017) Sahai, R., Vlemmings, W. H. T., Gledhill, T., et al. 2017, ApJ, 835, L13
  * Sánchez Contreras et al. (2019) Sánchez Contreras, C., Báez-Rubio, A., Alcolea, J., et al. 2019, A&A, 629, A136
  * Shu et al. (1994) Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781
  * Smartt (2009) Smartt, S. J. 2009, ARA&A, 47, 63.
  * Soker (1998) Soker, N. 1998, ApJ, 496, 833
  * Soker (1990) Soker, N. 1990, AJ, 99, 1869
  * Soker (2001) Soker, N. 2001, MNRAS, 325, 584
  * Soker (2002) Soker, N. 2002, ApJ, 568, 726
  * Soker (2020) Soker, N. 2020, Galaxies, 8, 26
  * Soker (2022) Soker, N. 2022, Research in Astronomy and Astrophysics, 22, 122003
  * Soker & Kashi (2012) Soker, N., & Kashi, A. 2012, ApJ, 746, 100
  * Soker & Kashi (2016) Soker, N., & Kashi, A. 2016, MNRAS, 462, 217
  * Soker & Livio (1994) Soker, N., & Livio, M. 1994, ApJ, 421, 219
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253
  * Swings (1974) Swings, J. P. 1974, A&A, 34, 333
  * Torres et al. (2018) Torres, A. F., Cidale, L. S., Kraus, M., et al. 2018, A&A, 612, A113
  * Vink (2015) Vink, J. S. 2015, Very Massive Stars in the Local Universe, 412, 77
  * Vink & Sander (2021) Vink, J. S., & Sander, A. A. C. 2021, MNRAS, 504, 2051
  * Zickgraf (2006) Zickgraf, F.-J. 2006, Stars with the B[e] Phenomenon, 355, 135
  * Zickgraf et al. (1996) Zickgraf, F.-J., Humphreys, R. M., Lamers, H. J. G. L. M., et al. 1996, A&A, 315, 510
  * Zickgraf et al. (1986) Zickgraf, F. J., Wolf, B., Stahl, O., Leitherer, C., & Appenzeller, I. 1986, A&A, 163, 119
  * Zou et al. (2020) Zou, Y., Frank, A., Chen, Z., et al. 2020, MNRAS, 497, 2855