11institutetext: Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer
de Can Magrans s/n, 08193 Barcelona, Spain
11email<EMAIL_ADDRESS>22institutetext: Ural Federal University, 19
Mira Str., 620002 Ekaterinburg, Russia 33institutetext: Research Institute of
Physics, Southern Federal University, Rostov-on-Don 344090, Russia

# The backreaction of stellar wobbling on accretion discs of massive
protostars

D. M.-A. Meyer 11 E. Vorobyov 2233

###### Abstract

Context.  In recent years, it has been demonstrated that massive stars see
their infant circumstellar medium shaped into a large, irradiated,
gravitationally unstable accretion disc during their early formation phase.
Such discs constitute the gas reservoir in which nascent high-mass stars gain
substantial fraction of their mass by episodic accretion of dense gaseous
circumstellar clumps, simultaneously undergoing accretion-driven bursts and
producing close-orbit spectroscopic companions to the young high-mass stellar
object.

Aims.  We aim to evaluate the effects of stellar motion, caused by the disc
non-axisymmetric gravitational field, on the disc evolution and its spatial
morphology. In particular, we analyze the disc propensity to gravitational
instability and fragmentation, and also disc appearance on synthetic
millimeter-band images pertinent to the alma facility.

Methods.  We employed three-dimensional radiation-hydrodynamical simulations
of the surroundings of a young massive star in the non-inertial spherical
coordinate system, adopting the highest spatial resolution to date and
including the indirect star-disc gravitational potential caused by the
asymmetries in the circumstellar disc. The resulting disc configurations were
postprocessed with the radiation transfer tool RADMC-3D and CASA softwear to
obtain disc synthetic images.

Results.  We confirm that the early evolution of the accretion disc is notably
different when stellar wobbling is taken into account. The redistribution of
angular momentum in the system makes the disc smaller and rounder, reduces the
number of circumstellar gaseous clumps formed via disc gravitational
fragmentation, and prevents the ejection of gaseous clumps from the disc. The
synthetic predictive images at millimeter wavelengths of the accretion disc
including stellar wobbling are in better agreement with the observations of
the surroundings of massive young stellar objects, namely,  AFGL 4176 mm1,
G17.64+0.16 and G353.273, than our numerical hydrodynamics simulations
omitting this physical mechanism.

Conclusions.  Our work confirms that stellar wobbling is an essential
ingredient to account for in numerical simulations of accretion discs of
massive protostars.

###### Key Words.:

giant planet formation – $\kappa$-mechanism – stability of gas spheres

##  1 Introduction

Although the formation of massive stars is a rare event, those objects are
preponderant engines in the cycle of life of matter in galaxies. The mystery
of their birth principally lies in the high optical depth of the parent pre-
stellar cores in which they form. Apart from revealing their common
multiplicity (Kraus et al., 2017), the quest for the first observation of
Keplerian-orbiting material around massive young stellar objects (MYSOs) saw a
first breakthrough with the works of Johnston et al. (2015); Ilee et al.
(2016); Forgan et al. (2016); Ginsburg et al. (2018); Maud et al. (2018). The
study of MYSOs continued with the hunt for evidence of filamentary/clumpy
substructures in discs around W33A MM1-Main (Maud et al., 2017) G350.69-0.49
(Chen et al., 2017) and G11.92-0.61 MM1 (Ilee et al., 2018). Besides,
accretion variability onto high-mass protostars has been monitored in massive
young stellar objects (Keto & Wood, 2006; Stecklum et al., 2017) and in pulsed
bipolar outflows (Cunningham et al., 2009; Cesaroni et al., 2010; Caratti o
Garatti et al., 2015; Purser et al., 2016; Reiter et al., 2017; Burns et al.,
2017; Burns, 2018; Purser et al., 2018; Samal et al., 2018).

Since the variability of the accretion flow onto the surface of massive
protostars governs their pre-main-sequence evolution in the Hertzsprung-
Russell diagram (Meyer et al., 2019a), the study of massive young stellar
objects is therefore tied to the understanding of the inhomogeneities and
substructures in their circumstellar discs as likely causes of the accretion
variability. These non-axisymmetric disc features are possibly produced by
instabilities of self-gravitating nature, which were extensively studied in
the context of low-mass stellar objects (Papaloizou & Savonije, 1991; Pickett
et al., 2003; Rafikov, 2005, 2007). Several questions remain wide open, such
as MYSOs chemical evolution (Ahmadi et al., 2018; Guadarrama et al., 2023) and
the role of magnetic fields in the regulation of their fragmentation
(Hennebelle et al., 2016).

The circumstellar medium of massive protostars has been investigated using
computationally-intensive numerical simulations, see e.g. Krumholz et al.
(2007a, b, 2009a); Harries (2015); Harries et al. (2017). Such calculations
raised a number of questions and revealed technical challenges caused by the
three-dimensional intrinsic nature of the problem and by the multiple physical
processes involved in the formation mechanism of massive stars (cooling and
heating of the gas, protostellar irradiation and radiation transport into the
disc, disc wind, non-ideal magneto-hydrodynamical processes, etc.). Amongst
many approaches developed to tackle the problem of massive star formation, the
new paradigm known as the burst mode of accretion produced results consistent
with not only models and observations of the lower mass regimes of star
formation, but also with direct observations of forming massive stars.

The burst mode of accretion is a depiction of star-forming processes which
fills the gap between the episodic accretion onto protostars and the luminous
bursts originating from young stellar objects and observed in some star-
forming regions (Vorobyov & Basu, 2005, 2006, 2010; Machida et al., 2011;
Vorobyov & Basu, 2015). In this picture, clumpy circumstellar material from
gravitationally-unstable accretion disc around young stars resulting from the
infall of the molecular material of pre-stellar cores inward-migrate from the
disc to the protostellar surface, hence triggering a sudden luminosity rise.
This provides an explanation to the so-called FU-Orionis accretion-driven
burst phenomenon and offers evident connections to the pre-zero-age-main-
sequence stellar evolution (Elbakyan et al., 2019) and the
prediction/observations of substructures in discs (Dong et al., 2016). The
general idea of the burst mode of accretion was extended and elaborated by
Nayakshin and co-authors to better explain the burst phenomenon and include
planet formation in inward-migrating clumps (Nayakshin, 2010; Nayakshin &
Lodato, 2012; Nayakshin, 2016, 2017; Vorobyov & Elbakyan, 2018; Elbakyan et
al., 2021, 2023).

Interestingly, this picture originally developed for the low-mass regime of
star formation can be naturally extended to the primordial and massive regimes
of star formation. In Meyer et al. (2017) and Meyer et al. (2018) it was shown
that a self-consistent picture unifying disc fragmentation and binary
formation in the context of massive young stellar objects is possible, linking
both accretion-driven outbursts as monitored in the MYSOs S255IR-NIRS3 and NGC
6334I-MM1 to the close/spectroscopic companions of some massive stars (Mahy et
al., 2013; Kobulnicky et al., 2014; Chini et al., 2012). Particularly,
accretion-driven outbursts are interpreted as to be the observational
signature of the presence of a self-gravitating disc shaped by efficient
gravitational instability. The bursts occur in series of multiple eruptions
and their intensity is a function of the mass of the accreted material (Meyer
et al., 2019c). Such episodic flares witness the sudden increment of the
protostellar mass with dense, compact circumstellar material (Elbakyan et al.,
2023) that affect the internal structure of the protostars, profoundly
changing their surface properties and provoking excursions in the Hertzsprung-
Russel diagram (Meyer et al., 2019a).

Table 1:  List of the simulation models performed in our study. The table provides the initial mass of the molecular pre-stellar core $M_{\rm c}$ (in $M_{\odot}$), its initial rotational-to-gravitational energy ratio $\beta$ (in $\%$) and the final simulation time $t_{\rm end}$ in each model. The last column indicates whether the simulation includes stellar wobbling or not.  $\mathrm{Models}$ |  $M_{\rm c}$ $(M_{\odot})$ | $\mathrm{Grid}\,\mathrm{resolution}$ |  $\beta$ ($\%$) |  $t_{\rm end}$ ($\rm kyr$) | $\mathrm{Wobbling}$
---|---|---|---|---|---
Run-512-100$\rm M_{\odot}$-4$\%$-wio | $100$ | $512\times 81\times 512$ | $4$ | 32.65 | no
Run-512-100$\rm M_{\odot}$-4$\%$-wi | $100$ | $512\times 81\times 512$ | $4$ | 32.65 | yes

Numerical results have preceded the direct observation of an accretion disc
around a young high-mass star. Such predictions divide into two distinct
categories of models, obtained either using a Cartesian grid with an adaptive
mesh refinement method and a protostar allowed to move freely inside the
computational domain (Krumholz et al., 2007a, 2009b; Peters et al., 2010;
Commerçon et al., 2011; Seifried et al., 2011, 2012, 2013; Klassen et al.,
2014; Seifried et al., 2015; Klassen et al., 2016; Rosen et al., 2016; Mignon-
Risse et al., 2021; Commerçon et al., 2022; Mignon-Risse et al., 2023a), or
using a spherical grid with a static logarithmic mesh and assuming that the
growing protostar is fixed to the origin of domain (Kuiper et al., 2010a, b,
2011; Kuiper & Yorke, 2013a, b; Meyer et al., 2017, 2018, 2019c, 2019b; Ahmadi
et al., 2019; Oliva & Kuiper, 2020; Meyer et al., 2021; Oliva & Kuiper, 2023a,
b). In both cases, the initial conditions and the included microphysical
processes in the models are similar (self-gravity, radiation transport coupled
to proto-stellar evolution).

Cartesian grids have the advantage not to have coordinate singularities and to
permit a straightforward implementation of the stellar motion, however, the
inclusion of stellar radiation transfer into the models makes calculations
computationally expensive. Spherical grids allow us to reach high spatial
resolution in the inner circumstellar environment of protostars, which
precisely permitted to reveal the burst mode of accretion in the context of
massive star formation (Meyer et al., 2017). Both computing philosophies have
been thoroughly compared on the basis of high-spatial resolution simulations
in Mignon-Risse et al. (2023b) and the authors conclude therein in a
qualitative agreement between both approaches, the remaining differences in
the results being minor.

Nevertheless, we showed in Meyer et al. (2022) that when permitting the
protostar to move, as a consequence of the disc-star gravitational
interaction, stellar motion delays disc fragmentation and affects the
millimeter continuum emission properties (Meyer et al., 2022). Such an effect
is not encompassed in the comparison study of Mignon-Risse et al. (2023b).
Difference between the spherical models of Meyer et al. (2022) and the
Cartesian simulations of Mignon-Risse et al. (2023b) lie in the spatial
resolution and the absence/presence of stellar wobbling. It is hereby proposed
to perform a comparison of disc simulation with the resolution of Mignon-Risse
et al. (2023b), with and without the inclusion of stellar inertia.

This study aims at investigating the effects of finite stellar inertia onto
the formation and growth of circumstellar accretion discs around massive
protostars. We also address the question of the millimeter-observability by
continuum dust emission of such disc nebulae and of the substructures forming
in them by gravitational instability. To this end, a pair of three-dimensional
gravito-radiation-hydrodynamics numerical simulations is performed, resolved
with the highest to-date spatial resolution used in spherically-based
calculation in the study of Mignon-Risse et al. (2023b). Two scenarios are
considered: without and with stellar inertia, namely, without and with the
star motion in response to the gravitational interaction with the disc.
Selected characteristic outcomes of the simulations are post-processed with
radiation transfer and imaging methods in order to produce synthetic $1.3\,\rm
mm$ dust continuum images for the Atacama Large Millimeter/submillimeter Array
(alma) interferometer operating in its Band 6, see the work of Meyer et al.
(2018) and Jankovic et al. (2019). We explore in these images the indirect
effects of stellar inertia onto the disc structure, its eventual fragmentation
by efficient gravitational instability, the formation of gaseous clumps
therein, as well as their dynamics inside (migration) and outside (ejection)
of the disc. Finally, the most realistic high-resolution disc model,
possessing both a high-spatial resolution and stellar inertia included in it,
is compared to observations available in the literature.

This work is organized as follows. In Section 2, the methods used to perform
high-resolution hydrodynamical simulations of a realistic accretion disc
around a young massive star are presented, together with the radiation
transfer procedure utilised to predict its continuum thermal dust emission at
the milliliter waveband. Section 3, introduce the reader to the results and
details of the hydrodynamical simulations of accretion discs surrounding
massive protostars. The radiative transfer calculations for dust continuum
emission of the stellar surroundings of these forming massive stars are shown
in Section 4, for several characteristic time instances of the disc evolution.
Conclusions are finally provided in Section 5.

##  2 Method

In this section, the methods used to simulate and predict the millimeter
continuum appearance as seen by the ALMA interferometre of the circumstellar
medium of a forming massive protostar are presented.

###  2.1 Radiation-hydrodynamics simulations

Numerical simulations of accretion discs around protostars were carried out
under the assumption of a midplane-symmetric computational domain. It is
initialised with a $M_{\rm c}=100\leavevmode\nobreak\ M_{\odot}$ pre-stellar
core of uniform temperature $T_{\rm c}=10\,\rm K$ which is rigidly rotating.
The initial mass density distribution is spherically symmetric and profiled
as,

$\rho(r)=\frac{3}{8\pi}\frac{M_{\rm c}}{R_{\rm c}^{3/2}}r^{-3/2},$ (1)

where $r$ is the radial coordinate and $R_{\rm c}$ the outer core radius,
respectively. The gas angular velocity of the pre-stellar core is initialised
according to the distribution $\Omega(R)\propto R^{-3/2}$, with the
cylindrical radius $R=r\sin(\theta)$, and according to the rotational-to-
gravitational energy ratio,

$\beta=\frac{E_{\rm rot}}{E_{\rm grav}}=0.04\propto\frac{1}{7}\frac{1}{GM_{\rm
c}R_{\rm c}^{-3}}\int_{0}^{\pi}d\theta\sin(\theta)^{3}.$ (2)

The grid mesh $[r_{\rm in},R_{\rm c}]\times[0,\pi/2]\times[0,2\pi]$ mapping
the computational domain expands logarithmically along the radial direction
$r$, as a cosine in the polar direction $\theta$ and is uniform along the
azimuthal direction $\phi$. It is made of $N_{\rm
r}=512\times\,N_{\rm\theta}=81\times\,N_{\rm\phi}=512$ grid zones,
respectively. The inner radius $r_{\rm in}$ constitutes a semi-permeable sink
cell fixed onto the origin of the domain and the outer radius, assigned to
outflow boundary conditions, is located at $R_{\rm c}=0.1\,\rm pc$. This study
adopts $r_{\rm in}=20\,\rm au$, which permits the highly-resolved models to
reach long integration times $t_{\rm end}$ without dealing with dramatic time-
step restrictions in the innermost grid zones. The utilised grid spatially
resolves the inner region of the midplane where the disc fragments are
accreted, while keeping the overall number of grid cells to a decent value and
thus reducing the numerical cost of this computationally-intensive
calculation.

The numerical simulations follow the gravitational collapse of the pre-stellar
core and the early evolution of circumstellar disc around the massive
protostar, accounting for both the stellar direct irradiation of protostellar
photospheric photons and the radiation transport into the circumstellar
medium. The accretion rate is calculated onto the protostar as the material
loss $\dot{M}$ through the sink cell and the surface and intrinsic properties
of the protostar, such as the stellar radius and the photospheric luminosity,
are time-dependently interpolated using the pre-calculated protostellar
evolutionary tracks of Hosokawa & Omukai (2009), see also Meyer et al. (2017,
2018). Two models are performed with the above described initial conditions,
with and without the inclusion of stellar inertia (Table 1).

###  2.2 Governing equations

The set of equations describing the dynamics of the pre-stellar core collapse
reads,

$\frac{\partial\rho}{\partial t}+\@vec{\nabla}\cdot(\rho\@vec{v})=0,$ (3)
$\frac{\partial\rho\@vec{v}}{\partial
t}+\@vec{\nabla}\cdot(\rho\@vec{v}\otimes\@vec{v})+\@vec{\nabla}p=\@vec{f},$
(4) $\frac{\partial E}{\partial
t}+\@vec{\nabla}\cdot((E+p)\@vec{v})=\@vec{v}\cdot\@vec{f},$ (5)

which are the relations for the conservation of the mass, the momentum and the
energy of the infalling material. In the above relations, the gas density
reads $\rho$, the thermal pressure $p=(\gamma-1)E_{\rm int}$, the gas velocity
$\@vec{v}$ and $\gamma=5/3$ is the adiabatic index. Hence,

$E=E_{\rm
int}+\rho\frac{\@vec{v}^{2}}{2}=\frac{p}{(\gamma-1)}+\rho\frac{\@vec{v}^{2}}{2},$
(6)

is the total energy of the gas. The right-hand term of the momentum and energy
conservation laws are the force density vector,

$\@vec{f}=-\rho\@vec{\nabla}\Phi_{\rm tot}-\lambda\@vec{\nabla}E_{\rm
R}-\@vec{\nabla}\cdot\Big{(}\frac{\@vec{F_{\star}}}{c}\Big{)}\@vec{e}_{\rm
r},$ (7)

with the flux limiter $\lambda=1/3$, the thermal radiation energy density
$E_{\rm R}$, the radial unit vector $\@vec{e}_{\rm r}$, the protostellar
radiation flux $\@vec{F_{\star}}$ and the speed of light $c$. The quantity
$\Phi_{\rm tot}$ is the total gravitational potential of the gas and the star,
including the indirect potential as described below.

Figure 1:  Mid-plane density fields (in $\rm g\,\rm cm^{-3}$ ) in the
accretion disc hydrodynamical simulations performed without (left panels) and
with (right panels) stellar inertia. The figures are shown for several time
instances, namely $19.15$ (a,b), $19.85$ (c,d) and $23.25\,\rm kyr$ (e,f). The
initial kinetic-to-gravitational energy ratio of the collapsing molecular
cloud is $\beta=4\%$.
Figure 2:  As Fig. 1, for the time instances $24.35$ (a,b), $27.55$ (c,d) and
$32.65\,\rm kyr$ (e,f).

###  2.3 Stellar inertia

This numerical setup uses the implementation for stellar motion presented in
Regály & Vorobyov (2017); Hirano et al. (2017); Meyer et al. (2022). It is
implemented as an additional indirect potential, $\Phi_{\mathrm{wobbling}}$
together with its associated force $\@vec{F}_{\mathrm{wobbling}}$. This
indirect force is implemented as the following acceleration,

$\@vec{g}^{{}^{\prime}}+\frac{\@vec{F}_{\mathrm{disc/\star}}}{M_{\star}}=0,$
(8)

where the disc-to-protostar gravitational force reads,

$\@vec{F}_{\mathrm{disc/\star}}=-GM_{\star}\int_{\mathrm{disc}}\frac{\rho(r)\,dV}{r^{2}}\@vec{e}_{\mathrm{r}},$
(9)

with the protostellar mass $M_{\star}$, the gravitational constant G, the
radial unit vector $\@vec{e}_{\mathrm{r}}$ and the mass in a volume element
$dM_{\mathrm{disc}}(r)=\rho(r)dV$. Hence, the total gravitational potential
reads $\Phi_{\rm tot}=\Phi_{\rm sg}+\Phi_{\mathrm{wobbling}}$, where
$\Phi_{\rm sg}$ is the input from the disc, the envelope, and the star. The
purpose of this study consists in exploring the effects of the stellar inertia
onto the dynamics and morphology of accretion discs at a spatial resolution
exceeding that of the previous study of this series, see Meyer et al. (2022).

###  2.4 Numerical methods

The equations of gravito-radiation-hydrodynamics are solved with the pluto
code111http://plutocode.ph.unito.it/ (Mignone et al., 2007, 2012). The direct
irradiation feedback of the protostar and the radiation transport in the
accretion disc are both taken into account within the gray approximation. It
uses the numerical scheme presented in Kolb et al.
(2013)222http://www.tat.physik.uni-tuebingen.de/ pluto/pluto_radiation/ and
adapted in Meyer et al. (2018) for the study of massive star formation. It is
an algorithm which first ray-traces photon packages from the protostellar
surface to the disc and then diffuses their propagation into it in the flux-
limited approximation. This approach allows us to accurately consider both the
inner heating and the outer cooling of our irradiated accretion discs, see
Vaidya et al. (2011). We note that similar radiation-hydrodynamics methods are
also presented in Commerçon et al. (2011), Flock et al. (2013) and Bitsch et
al. (2014).

The opacity description as well as the estimate of the local dust properties
are similar as in Meyer et al. (2018), where the gas $T_{\rm gas}$ and dust
$T_{\rm dust}$ temperatures are calculated assuming the equilibrium between
the silicate grains temperature and the total radiation field. Stellar gravity
is modelled by calculating the total gravitational potential of the central
protostar and include the self-gravity of the gas by solving the Poisson
equation using the PETSC library333https://www.mcs.anl.gov/petsc/. We neglect
turbulent viscosity by assuming that the most efficient mechanism for the
transport of angular momentum are the gravitational torques in a self-
gravitating disc. The effects of the stellar inertia onto the disc wobbling in
included as described in Meyer et al. (2019c). The reader interested in
further reading about the method can refer to the other papers in this series,
initiated with Meyer et al. (2017).

###  2.5 Construction of synthetic disc images

As a diagnostic to evaluate the effects of stellar inertia on realistic, high-
resolution models of the circumstellar medium of massive protostars, synthetic
alma images of the accretion discs were calculated using the post-processing
method presented in Meyer et al. (2019b). For each panel displayed in Figs.
1-2 the dust density fields of the embedded accretion discs simulated with the
pluto code were imported into the radiative transfer code
radmc-3d444http://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/
(Dullemond, 2012). The standard dust to gas mass ratio of 0.01 was assumed
when converting the gas density into that of dust. Then, the dust temperature
is derived by Monte-Carlo calculation on the basis of the dust density, using
the method presented in Bjorkman & Wood (2001) and using $10^{10}$ photons
packages ray-traced from the protostellar surface to the outer region of the
accretion disc. The proper radiative transfer calculation against dust opacity
of the star light scattered into the dusty disc is subsequently performed for
a monochromatic light of wavelength centered onto $1.2\,\rm mm$ ($249.827\,\rm
Ghz$ with a channel width of $50.0\,\rm Mhz$), considering that the dust is a
Laor & Draine (1993) mixture of silicate particles. The photopshere is
modelled as a black body of effective temperature $T_{\rm eff}$ that is
derived for each selected simulation snapshots as a function of the age of the
protostar, i.e. its mass and current accretion rate, which are used to
bilinearly interpolate the protostellar evolutionary tracks of Hosokawa et al.
(2010).

Finally, synthetic images of the accretion discs were generated using the
Common Astronomy Software Applications casa555https://casa.nrao.edu/ (McMullin
et al., 2007) using the radiation flux output from radmc-3d. No inclination
angle for the accretion disc with respect to the plane of the sky is assumed
throughout the whole simulation process and images were produced with a field
of view of $2000\,\rm au$ around the growing protostar. The simulated alma
interferometric images were thus obtained. They can either serve as a test for
the observability of the simulated accretion discs, or as a tool to directly
compare the disc models with available and/or forthcoming observations to be
performed by the alma facility.

As in Meyer et al. (2019b), it is assumed that the conditions for the
synthetic observations are ideal. In other words, the pattern of
telescopes/antennae are considered to be in their most extended spatial
configuration at the Llano de Chajnantor plateau, permitting the acquisition
of long-baseline images resolved with the maximal possible spatial resolution
with the smallest beam size. This resolution is about $0.015^{"}$ for the
C43-10 configuration that makes use of 43 $12\,\rm m$ antennae. The modelling
of these Cycle 10 synthetic alma observations assumes an exposure time of $10$
minutes and a precipitable water vapour $\rm pwv=0.6$. The assumed distance to
the source is taken to $1\,\rm kpc$, which corresponds to that of the Earth’s
closest high-mass star-forming region Orion.

Figure 3:  Synthetic $1.2\,\rm mm$ dust continuum emission maps as seen by the
alma  interferometer in its antenna configuration 10, of the accretion disc
simulated without (left panels) and with (right panels) stellar inertia. The
panels are shown for several time instances displayed in Fig. 1, namely,
$19.15$ kyr (a,b), $19.85$ kyr (c,d) and $23.25\,\rm kyr$ (e,f). The distance
to the source is assumed to be $1\rm kpc$.
Figure 4:  As Fig. 3, for the time instances displayed in displayed in Fig. 2,
namely $24.35$ kyr (a,b), $27.55$ kyr (c,d) and $32.65\,\rm kyr$ (e,f).

##  3 Results

This section presents the performed high-resolution disc models and details
their accretion disc morphological properties.

###  3.1 Disc evolution: $19.15-23.25\,\rm kyr$ after onset of the collapse

Fig.1 displays the midplane density field of the accretion disc (in $\rm
g\,\rm cm^{-3}$) in our simulations $\rm Run$$-$$512$$-$$100\rm
M_{\odot}$$-$$4$%$-\rm wio$ (without stellar wobbling, left-hand panels) and
$\rm Run$$-$$512$$-$$100\rm M_{\odot}$$-$$4$%$-\rm wi$ (with stellar wobbling,
right-hand panels). The discs models are plotted at several characteristic
time instances such as $19.15\,\rm kyr$ (top series of panels), $19.85\,\rm
kyr$ (middle series of panels), $23.25\,\rm kyr$ (bottom series of panels)
after onset of the gravitational collapse of the molecular cloud. The panels
depict the inner ($\leq 1600\,\rm au$) region of the computational domain. The
simulation begin with the gravitational collapse of the rotating pre-stellar
molecular cloud and an accretion disc forms when the infalling flow hits the
centrifugal barrier near inner computational boundary engendered by the
conservation of angular momentum of the infalling matter.

At time $19.15\,\rm kyr$, accretion discs in both models have a round
morphology with radii of about a few $100\,\rm au$. The accretion disc is
already structured with several spiral arms wrapped around the inner region of
the disc. At that early time, differences between the models are noticeable in
the sense that the disc modelled without stellar inertia is larger by almost a
factor $\sim 2$ and the spiral arms are more pronounced, extending to larger
radii (Fig.1a), while the disc in the model with stellar inertia is still
smaller (Fig.1b). A gaseous clump has formed by efficient gravitational
instability in the disc without stellar inertia, which is the sign of faster
gravitational fragmentation of the disc when the central protostar is not
allowed to move (Fig.1a).

At time $19.85\,\rm kyr$, one can see that the disc modelled without stellar
inertia keeps on growing in size up to a radius $\sim 450\,\rm au$. Its
internal structure includes spiral arms of different sizes and thicknesses,
separated by a more diluted interarmed space. The innermost arm is thinner,
denser and hosts a circumstellar gaseous clump that is currently migrating
downward to the inner disc region, where the protostar is located and accretes
mass. This mechanism is responsible for the production of accretion-driven
outbursts (Meyer et al., 2017) and/or the formation of close binary companions
to the young high-mass star (Meyer et al., 2018). Another clump, also formed
by gravitational instability in the accretion disc, is ejected by the
gravitational swing three-body mechanism away from its parent accretion disc
(Fig.1c). This is a major difference between the simulations without stellar
inertia and those in which stellar inertia is included. Indeed, the disc
modelled with stellar wobbling continues to display its round, slightly
structured morphology, but it does not exhibit any sign of fragmentation up to
this time instance (Fig.1d). This confirms, on the basis of models with higher
spatial resolution, the delay of disc fragmentation induced by stellar inertia
compared to simulations neglecting the stellar motion in response to the disc-
star gravitational interaction (Meyer et al., 2022).

At $23.25\,\rm kyr$, the model $\rm Run$$-$$512$$-$$100\rm
M_{\odot}$$-$$4$%$-\rm wio$ that does not include stellar inertia displays an
accretion disc, which has violently fragmented and exhibits all associated
features, such as several spiral arms, many gaseous circumstellar clumps
(Meyer et al., 2018), as well as a trail left behind the ejected clump (see
Fig.1c). We note that this clump has now travelled away from the disc at a
distance $\geq 1600\,\rm au$ (Fig.1e). The clump ejection mechanism was
described for the first time in the context of low-mass star formation, when
gravitationally fragmenting accretion discs of solar-type stars generate free-
floating clumps thrown through the ISM, see the studies of Basu & Vorobyov
(2012); Vorobyov (2016). This time instance also reveals the first clear sign
of efficient gravitational instability in the accretion disc model simulated
with stellar inertia. A characteristic disc morphology develops, made of
several spiral arms of similar extent reaching $\geq 800\,\rm au$, together
with the establishment of a dense disc region where both spiral arms originate
and gaseous clumps migrate (Fig.1f). The two numerical models differ only
slightly by the disc radius at this time instance, but significantly by the
number of gaseous clumps in them, which is notably higher in the case ignoring
stellar inertia. This strengthens the conclusions of Meyer et al. (2022)
regarding the slow down and protracting effects of disc wobbling in the
development of gravitational instability in accretion discs around massive
protostars (Fig.1e).

###  3.2 Disc evolution: $24.35-32.65\,\rm kyr$ after onset of the collapse

Fig.2 presents the disc evolution at later times. In particular, the gas
spatial distributions are plotted at times $24.35\,\rm kyr$ (top series of
panels), $27.55\,\rm kyr$ (middle series of panels), $32.65\,\rm kyr$ (bottom
series of panels) after onset of the gravitational collapse of the molecular
cloud. The disc without stellar wobbling continues fragmenting and becoming
larger as it accretes mass, extending its denser spiral arms and forming
several disc-like structures around the clumps that are located at large
orbits ($\sim 800\,\rm au$), see Fig.2a. Most important differences between
the two models are the disc radial extend, as well as the level of
fragmentation of the disc. In particular, the disc with stellar inertia
included is more compact and has smaller clumps not surrounded yet by their
own accretion discs, see Fig.2b. As underlined in Meyer et al. (2022) and in
Regály & Vorobyov (2017) in the context of low-mass stars, the inclusion of
stellar inertia in the simulations permit the transfer of part of the disc
angular momentum to the wobbling star. Consequently, the disc angular momentum
is smaller compared to the case without stellar inertia. Since disc
fragmentation takes place at larger radii, compacter discs are expected to
fragment less (Johnson & Gammie, 2003).

Later, at times $24.35-27.55\,\rm kyr$ the differences between accretion disc
modelled with and without stellar inertia augment and become more pronounced.
The disc with motionless protostar keeps on expanding up to a radius $\sim
1500\,\rm au$, and most gaseous clumps themselves are surrounded by dense
disc-like structures (Fig.2c). Those structures around clumps in their turn
experience accretion and acquire mass from the disc material, while their
cores become denser and hotter, potentially up to the temperature of molecular
hydrogen dissociation, see Meyer et al. (2018, 2022). A smaller size of the
disc modelled with stellar inertia weakens the development of gravitational
fragmentation, and a sole migrating clump on the way to be accreted by the
protostar is visible in the inner disc region (Fig.2d).

Finally, at time $32.65\,\rm kyr$, the disc with the fixed central young high-
mass star has the most extended and complex morphology, composed of a pattern
of spiral arms of different densities and lengths, interspersed with gaseous
clumps, themselves potentially surrounded by nascent accretion discs (Fig.2e).
This picture is a scaled-up version of the development of massive
protoplanetary discs around young solar-type stars, see Vorobyov & Basu
(2010), and in the context of primordial star formation Vorobyov et al.
(2013). At this time instance, the accretion disc model with stellar inertia
included exhibits clear signs of fragmentation, with a large-scale spiral arm
enrolled around the disc structure, and several dense clumps migrate towards
the inner circumstellar region (Fig.2f). The disc size is similar to that of
the disc with fixed protostar when it started fragmenting violently (Fig.2a).
We note the still-ongoing displacement of the massive protostar away from its
initial location at the beginning of the gravitational collapse, when it
coincided with the geometrical centre of the pre-stellar core, at the origin
of the computational domain.

##  4 Observational manifestations of stellar wobbling

This section presents predictions for the appearance of the accretion discs of
young massive stars that we simulated, as seen by the ALMA interferometer. It
compares the synthetic images and draws conclusions regarding the necessary
and relevant physics to be included in models when simulating discs around
massive protostars.

###  4.1 Synthetic ALMA images of the accretion discs

Fig.3 shows the synthetic millimetre images of the accretion discs as seen by
the alma facility. The plots distinguish models without (left-hand side
panels) and with (right-hand side panels) and the discs models are displayed
at several characteristic time instances, namely $19.15\,\rm kyr$ (top series
of panels), $19.85\,\rm kyr$ (middle series of panels), $23.25\,\rm kyr$
(bottom series of panels) after onset of the gravitational collapse of the
molecular cloud. The panels display the inner ($\leq 1600\,\rm au$) region of
the computational domain, which, at a distance of $1\,\rm kpc$ corresponds to
a region in the plane of the sky that is $1.0\,\rm arcsec$ wide, in both
longitude and latitude, respectively. Both disc synthetic images at time
$19.15\,\rm kyr$ have the appearance of a bright spot, representing the
nascent accretion disc which has just formed right after the end of the free-
fall gravitational collapse (Fig.3a,b). Nothing, except for two tiny
protuberances, giving the image of the disc modelled with stellar inertia a
slightly eared-like structure (Fig.3a), distinguish the two models (Fig.3b).
We note that the images that we produce do not assume any inclination of the
disc regarding to the plane of the sky. Differences arise at time $19.85\,\rm
kyr$, when the disc without stellar inertia initiates fragmentation and
develops a large spiral arm, which is enrolled all around the star and has a
clump in it. Both substructures are observable by alma (Fig.3c), whereas the
image of the disc with stellar inertia, which is yet not fragmented, conserves
its round morphology (Fig.3d).

At time $23.25\,\rm kyr$, all accretion disc models have undergone
gravitational fragmentation, however, their thermal dust emission at $1.3\,\rm
mm$ is quantitatively similar, despite of clear morphological differences of
the disc density fields. In both models, the infrared disc nebula reveals a
bright central region of global circular shape, surrounded by bright spots
corresponding to the migrating gaseous clumps (Fig.3e,f). In the case of the
observation of such discs, both will conclude on the presence of a fragmenting
disc around a massive protostar. It is interesting to note that the flying
clump catapulted from the accretion disc modelled without stellar inertia
(Fig.1c) is not visible in the infrared, as the synthetic image of Fig.3c,e
demonstrates. This can be explained by the poorer numerical resolution of the
outer regions of the computational domain compared to the inner disc region.
Such a lack of resolution has long-time been the reason why local numerical
simulations of massive star formation could not consistently conclude on disc
fragmentation (Krumholz et al., 2009b). The density and temperature contrast
across the ejected clump washes out as it enters the outer regions with
progressively lower numerical resolution, making it more difficult to detect.
This possible explanation is strengthened by the images at time $24.35\,\rm
kyr$, still displaying two northern bright spots corresponding to two massive
clumps (hence dense and hot) surrounded by accretion structures in the disc
(Fig.2a), while the ejected lower-mass clump is not visible in the field-of-
view.

Fig.4 plots the later evolution of the disc as seen by ALMA at times
$24.35\,\rm kyr$ (top series of panels), $27.55\,\rm kyr$ (middle series of
panels), $32.65\,\rm kyr$ (bottom series of panels). It essentially indicates
that the increasing complexity of the accretion disc simulated without stellar
inertia does not translate entirely into the millimeter images. For example,
the very complex circumstellar patterns of Fig.1, 2 are not fully recognisable
in Fig.4c but also in Fig.4e. Nevertheless, we note that a circumstellar disc
around a gaseous clump located in the southern part of the disc in Fig.4c
appears as a circular, enrolled structure in Fig.4c. This indicates that the
surroundings of circumstellar clumps that are in their turn accreting disc
material, i.e. are on the way to low-mass star formation and can be considered
as embryos of nascent secondary stars, could be observed with today’s
available facilities. The disc substructure are barely visible at later time
$32.65\,\rm kyr$ in the model without stellar inertia, except in the vicinity
of the protostar (Fig.4e). On the other hand, the synthetic image of the
accretion disc modelled with stellar inertia is well traced by the emission,
since it is denser, and therefore it has accumulated more dust mass from the
infalling material which reemits the protostellar light more efficiently and
the discs are easier to be observed by alma.

###  4.2 Caveats of simulation method

The limitations of the methods used to produce the models presented in this
study remain essentially the absence of the pre-stellar core magnetic field in
the initial conditions, directly associated to the difficulty to include non-
ideal effects, such as Ohmic diffusion and Hall effect, in the early disc
formation, as well as the use of a radiation transport scheme restricted to
non-ionizing photons. No disc wind or photoevaporation mechanism are also
present in the model. These caveats remain sensibly similar to those already
presented in the precedent papers of this series, see for example Meyer et al.
(2018). The most important improvement would be to include the magnetic field
in the disc simulations, since the photoionizing feedback is principally
relevant for the physics of the bipolar region growing perpendicularly to the
circumstellar disc. Only then, the construction of sophisticated simulated
disc images fine-tuned to particular objects will be possible. Apart from the
physical processes that are absent in our numerical model, potential future
improvements could be the use of a computational grid without the mid-plane
symmetry that is currently utilised in the present study and the adoption of a
sink cell radius of value $<20\,\rm au$. This will better capture any vertical
oscillations of the accretion disc and study the vertical effects of non-
aligned disc of massive protostars. Such update would nonetheless solely be
possible at the cost of a reduced time-step and therefore of a dramatic
increase of the computational costs of the numerical simulations.

###  4.3 Comparison with previous numerical works

The work of Mignon-Risse et al. (2023b) performs a comparison between the
ramses and pluto codes regarding the initial gravitational collapse of a
massive pre-stellar core, the first fragmentation era, and the disk
fragmentation of the circumstellar disc based on the initial conditions and
grid discretization of the computational domain corresponding to the highest
resolution simulation models of Oliva & Kuiper (2020). Both codes include
similar microphysical processes, such as radiation transport with direct
stellar irradiation, central gravity, and self-gravity of the gas. Sub-grid
modules are included for disc fragmentation via sink particles, stellar
evolution, and evolution of the dust component of the circumstellar material.
The spatial resolution is the highest ever achieved for calculating disc
fragmentation around a nascent massive star (Oliva & Kuiper, 2020). Very good
agreements between the two codes have been found during the early free-fall
collapse of the pre-stellar core and the era of the first fragmentation. An
overall satisfying match is reported in the properties of the disc fragments,
testifying to the consistency in the radiation-hydrodynamical context with
self-gravity of both ramses (Teyssier, 2002; Fromang et al., 2006) and pluto
(Mignone et al., 2007, 2012), despite differences in the numerical schemes
used. Indeed, pluto utilizes a ray-tracing method for direct irradiation with
disc radiation transport by Flux-Limited-Diffusion (Kolb et al., 2013), while
ramses is equipped with an M1 method (Rosdahl & Teyssier, 2015) for
protostellar feedback coupled to a Flux-Limited-Diffusion module (Mignon-Risse
et al., 2020) for disc radiation physics. ramses operates on a Cartesian grid
with adaptive mesh refinement, while pluto works with a spherical grid with
logarithmically-expanding spacing in the radial direction. We refer the
interested reader to further details on the algorithm differences as well as a
deep discussion on the advantages and conveniences of the calculation grids in
Mignon-Risse et al. (2023b).

The present paper compares two models computed with the pluto code, with an
implementation of the microphysical processes described in Meyer et al.
(2022). A significant difference is found between the models without and with
stellar inertia, which was previously reported in Meyer et al. (2022) based on
lower-resolution simulations. Disc wobbling acts as an additional effect
resulting in slowing down and prolonging the development of gravitational
instability in the accretion disc, reducing the number and magnitude of the
accretion-driven bursts. This is consistent with the comparison of the
accretion rate histories in Fig. 7 of Mignon-Risse et al. (2023b), which
reports, despite a similar quiescent accretion rate, much milder accretion
peaks onto the growing young massive star in the early $10\,\rm kyr$ of its
existence and indicates a smaller mass of the circumstellar clumps forming in
the disc. The difference in the dynamics of the disc fragments is interpreted
as originating from discrepancies in the central star mass due to the stellar
accretion models utilized in the two setups. Our results confirm the
morphological difference and changing fragment patterns in the disc as a
result of stellar wobbling, as noted by Meyer et al. (2022), with a spatial
resolution that matches that of both the ramses and pluto models in Mignon-
Risse et al. (2023b). We conclude that our results are consistent with the
discrepancies noted in Mignon-Risse et al. (2023b) and propose that stellar
wobbling is the factor explaining the differences between ramses and pluto
results.

###  4.4 Comparison with observations

The following section compares the modelled accretion discs with several
available observations of Keplerian structures around young massive stars.

####  4.4.1 The disc of AFGL 4176 mm1

The accretion disc surrounding the massive protostar AFGL 4176, also known as
G308.918+0.123 or IRAS 13395-6153, as documented by Persi et al. (1986),
exhibits circumstellar characteristics (Bøgelund et al., 2019) and follows a
Keplerian pattern of rotation (Johnston et al., 2020a). With a radius of
approximately $1000\,\rm au$, the disc features an asymmetric spiral arm
detected through alma 1.2 mm emission (Johnston et al., 2020b). The protostar
is $\approx 25\,\rm M_{\odot}$, its surrounding disc is $\approx 12.5\,\rm
M_{\odot}$ and it is located at a distance of $\approx 4.2\,\rm kpc$, see
Johnston et al. (2020a, b). The disc spiral arm extends up to about $1000\,\rm
au$ from the protostar and it has an inner ($\leq 900\,\rm au$) substructure
suggesting probable active disc fragmentation at work in the inner part. In
our simulations, structures of size reaching radii of $1000\,\rm au$ are
obtained at time $\geq 24.35\,\rm kyr$ when disc wobbling is excluded (Fig.
2(c) and (e)) and when gravitational fragmentation is very active in the disc.
The disc in the model with stellar wobbling is in general more compact and
features a spiral of comparable size to that of AFGL 4176 only at $t=32.65$
kyr (Fig. 2(f)). However, these differences become less evident when synthetic
images are compared (Fig. 4) and hence the model with wobbling can also be
consistent with the spiral features of AFGL 4176.

####  4.4.2 The disc of G17.64+0.16

The surroundings of the massive ($45\pm 10\,\rm M_{\odot}$) O-type protostar
G17.64+0.16 has been shown via alma long-baseline observations at band 6 of
$\rm H_{2}O$ emission to be made of a disc structure, which material is
characterized by Keplerian rotation around the central young star (Maud et
al., 2019). This disc reveals a strong concentration of dust particles
distributed in a ring-like outer region, possibly trapped by the pressure
field as transitional discs often do in the context of low-mass protostars, as
well as clear dense substructures in the inner disc region. The Toomre
analysis of the $\approx 100\,\rm pc$-large disc in G17.64+0.16 shows that no
gravitational instabilities are at work, and, hence, it is not in a state of
violent fragmentation. This apparent stability, despite the presence of
substructures, is puzzling and somehow challenges our results. The size of our
accretion disc without stellar wobbling exceeds that of G17.64+0.16, when the
first clump forms in the disc. On the other hand, such a large disc without
strong fragmentation at work is found in our model with stellar inertia (see
Fig. 1a,d). This makes stellar wobbling a necessary ingredient to reproduce
the qualitative features of the circumstellar medium of G17.64+0.16.

The particularly high mass of G17.64+0.16 is also not in accordance with our
series of models, see for example fig. 2 in the study of Meyer et al. (2022)
in which the mass of the protostar does not grow above $\approx 25\,\rm
M_{\odot}$. Since (i) a higher initial ratio of rotational-to-gravitational
energy results in the formation of protostars of lower mass but discs of
larger size, as demonstrated in the parameter study of Meyer et al. (2021),
and that (ii) the disc of G17.64+0.16 is rather small and young, we conclude
that the initial mass distribution of the G17.64+0.16 pre-stellar molecular
core nay have been drastically different from that in our models. Fig. 2g,h of
Meyer et al. (2021) shows that only molecular cores of initial mass $>180\,\rm
M_{\odot}$ can produce protostars entering the high-mass regime during the
core phase of free-fall gravitational collapse, while reaching $\approx
35\,\rm M_{\odot}$ in its direct follow-up evolutionary phase. We propose that
G17.64+0.16 has formed via the collapse of a massive
($>180\leavevmode\nobreak\ M_{\odot}$) but slowly rotating core.

####  4.4.3 The disc of G353.2731

High-resolution images with the Jansky-Very Large Array (J-VLA) and the
Australia Telescope Compact Array (ATCA) of the region of G353.2731 at radio
continuum and maser emission has shown that its circumstellar medium is shaped
as an accretion disc that is faced-on regarding to our observing line-of-
sight, see the study of Motogi et al. (2017). The mass of the protostar is
constrained to be about $10\,\rm M_{\odot}$, and it is surrounded by an
envelope. The inner region of the envelope, with a radius spanning distances
of approximately $\approx 100\,\rm au$, was interpreted in Motogi et al.
(2017) to be smaller than the typical disc sizes predicted by previous
numerical simulations. Indeed, the high-resolution models including radiative
transport and detailed stellar irradiation of Meyer et al. (2018) show that
spiral arms of such discs can easily reach $\approx 1000\,\rm au$ once
fragmentation is at work. This is confirmed by our high-resolution model
without stellar wobbling (Fig. 2c,e). Since the main effect of stellar inertia
that we highlight in our study is to delay disc fragmentation and to produce
discs of reduced size, we conclude that the small disc of G353.2731 is
consistent with our work.

##  5 Conclusions

This paper explores the effects of stellar inertia on the early development
and fragmentation of circumstellar accretion discs around a massive protostar.
Three-dimensional high-resolution numerical gravito-radiation-hydrodynamical
simulations are performed with the pluto code (Mignone et al., 2007, 2012;
Vaidya et al., 2018). The models begin at the onset of the gravitational
collapse of a rotating molecular cloud of kinetic-to-gravitational energy
ratio $\beta=4\%$ and are followed up to a time about $32\,\rm kyr$. The disc
models present the following two characteristic features: a high spatial
resolution corresponding to the most highly resolved simulation using a
spherical coordinate system (Mignon-Risse et al., 2023b) and a careful
treatment of the gravitational interaction at work between the disc and the
central young high-mass star (Meyer et al., 2022). The latter is included as
the acceleration the protostar undergoes when interacting gravitationally with
the non-axisymmteric accretion disc (Michael & Durisen, 2010; Hosokawa et al.,
2016; Regály & Vorobyov, 2017). The effects on the disc morphology of the
resulting stellar wobbling with respect to the coordinate centre is evaluated
by producing synthetic millimeter alma observations of the star-disc system,
which are discussed in the context of available observations.

Notable differences are found between the disc structure with and without
stellar wobbling. When the star is allowed to move, the development of
gravitational instabilities is slowed down and delayed. The disc keeps
displaying a small, round morphology, in which the formation of spiral arms
and gaseous clumps happens at later times compared to simulations performed
without wobbling. This confirms the findings of Meyer et al. (2022) obtained
on the basis of simulations with lower-resolution.  High-resolution accretion
discs modelled with stellar inertia generate a much simpler pattern. It
includes less gaseous clumps in it, around which secondary disc-like nebulae
do not form as easily and numerously as in the model with fixed central star.
Gaseous clumps ejected from the disc in the model without stellar wobbling
Basu & Vorobyov (2012); Vorobyov (2016) are not found once stellar inertia is
included into the model, at least for the kinetic-to-gravitational energy
ratio that we consider ($\beta=4\%$).

Synthetic $1.2\,\rm mm$ images of thermal dust emission of circumstellar discs
wherein stellar wobbling is considered turn out to be much more qualitatively
consistent with real observations of the discs of massive protostars than the
images produced in models without stellar inertia, in which the star is fixed
to the coordinate origin of the domain.  This applies to the observations of
the Keplerian circumstellar medium of AFGL 4176 mm1, G17.64+0.16 and G353.273.
This study stresses the importance of spatial resolution in the study of
accretion discs around massive protostars and it demonstrates that including
the star-to-disc gravitational interaction is a preponderant ingredient to
account for in realistic modelling of the surroundings of young massive
stellar objects.

###### Acknowledgements.

We are thankful to the anonymous referee for comments that helped to improve
the manuscript. The author acknowledge the North-German Supercomputing
Alliance (HLRN) for providing HPC resources that have contributed to the
research results reported in this paper. This research made use of the pluto
code developed at the University of Torino by A. Mignone
(http://plutocode.ph.unito.it/) and of the radmc-3d code developed at the
University of Heidelberg by C. Dullemond (https://www.ita.uni-
heidelberg.de/$\sim$dullemond/software/radmc-3d/). The figures have been
produced using the Matplotlib plotting library for the Python programming
language (https://matplotlib.org/). The data underlying this article will be
shared on reasonable request to the corresponding author. This work has been
supported by the grant PID2021-124581OB-I00 funded by
MCIN/AEI/10.13039/501100011033 and 2021SGR00426 of the Generalitat de
Catalunya. This work was also supported by the Spanish program Unidad de
Excelencia María de Maeztu CEX2020-001058-M. This work also supported by MCIN
with funding from European Union NextGeneration EU (PRTR-C17.I1). E.I.V.
acknowledges support from the Russian Science Foundation, project No.
23-12-00258.

## References

  * Ahmadi et al. (2018) Ahmadi, A., Beuther, H., Mottram, J. C., et al. 2018, A&A, 618, A46
  * Ahmadi et al. (2019) Ahmadi, A., Kuiper, R., & Beuther, H. 2019, A&A, 632, A50
  * Basu & Vorobyov (2012) Basu, S. & Vorobyov, E. I. 2012, ApJ, 750, 30
  * Bitsch et al. (2014) Bitsch, B., Morbidelli, A., Lega, E., & Crida, A. 2014, A&A, 564, A135
  * Bjorkman & Wood (2001) Bjorkman, J. E. & Wood, K. 2001, ApJ, 554, 615
  * Bøgelund et al. (2019) Bøgelund, E. G., Barr, A. G., Taquet, V., et al. 2019, A&A, 628, A2
  * Burns (2018) Burns, R. A. 2018, in IAU Symposium, Vol. 336, Astrophysical Masers: Unlocking the Mysteries of the Universe, ed. A. Tarchi, M. J. Reid, & P. Castangia, 263–266
  * Burns et al. (2017) Burns, R. A., Handa, T., Imai, H., et al. 2017, MNRAS, 467, 2367
  * Caratti o Garatti et al. (2015) Caratti o Garatti, A., Stecklum, B., Linz, H., Garcia Lopez, R., & Sanna, A. 2015, A&A, 573, A82
  * Cesaroni et al. (2010) Cesaroni, R., Hofner, P., Araya, E., & Kurtz, S. 2010, A&A, 509, A50
  * Chen et al. (2017) Chen, X., Ren, Z., Zhang, Q., Shen, Z., & Qiu, K. 2017, ApJ, 835, 227
  * Chini et al. (2012) Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925
  * Commerçon et al. (2011) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35
  * Commerçon et al. (2022) Commerçon, B., González, M., Mignon-Risse, R., Hennebelle, P., & Vaytet, N. 2022, A&A, 658, A52
  * Commerçon et al. (2011) Commerçon, B., Hennebelle, P., & Henning, T. 2011, ApJ, 742, L9
  * Cunningham et al. (2009) Cunningham, N. J., Moeckel, N., & Bally, J. 2009, ApJ, 692, 943
  * Dong et al. (2016) Dong, R., Vorobyov, E., Pavlyuchenkov, Y., Chiang, E., & Liu, H. B. 2016, ApJ, 823, 141
  * Dullemond (2012) Dullemond, C. P. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library
  * Elbakyan et al. (2023) Elbakyan, V. G., Nayakshin, S., Meyer, D. M. A., & Vorobyov, E. I. 2023, MNRAS, 518, 791
  * Elbakyan et al. (2021) Elbakyan, V. G., Nayakshin, S., Vorobyov, E. I., Caratti o Garatti, A., & Eislöffel, J. 2021, A&A, 651, L3
  * Elbakyan et al. (2019) Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146
  * Flock et al. (2013) Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43
  * Forgan et al. (2016) Forgan, D. H., Ilee, J. D., Cyganowski, C. J., Brogan, C. L., & Hunter, T. R. 2016, MNRAS, 463, 957
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Ginsburg et al. (2018) Ginsburg, A., Bally, J., Goddi, C., Plambeck, R., & Wright, M. 2018, ApJ, 860, 119
  * Guadarrama et al. (2023) Guadarrama, R., Vorobyov, E., Rab, C., et al. 2023, arXiv e-prints, arXiv:2312.07184
  * Harries (2015) Harries, T. J. 2015, MNRAS, 448, 3156
  * Harries et al. (2017) Harries, T. J., Douglas, T. A., & Ali, A. 2017, MNRAS, 471, 4111
  * Hennebelle et al. (2016) Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8
  * Hirano et al. (2017) Hirano, S., Hosokawa, T., Yoshida, N., & Kuiper, R. 2017, Science, 357, 1375
  * Hosokawa et al. (2016) Hosokawa, T., Hirano, S., Kuiper, R., et al. 2016, ApJ, 824, 119
  * Hosokawa & Omukai (2009) Hosokawa, T. & Omukai, K. 2009, ApJ, 691, 823
  * Hosokawa et al. (2010) Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478
  * Ilee et al. (2018) Ilee, J. D., Cyganowski, C. J., Brogan, C. L., et al. 2018, ApJ, 869, L24
  * Ilee et al. (2016) Ilee, J. D., Cyganowski, C. J., Nazari, P., et al. 2016, MNRAS, 462, 4386
  * Jankovic et al. (2019) Jankovic, M. R., Haworth, T. J., Ilee, J. D., et al. 2019, MNRAS, 482, 4673
  * Johnson & Gammie (2003) Johnson, B. M. & Gammie, C. F. 2003, ApJ, 597, 131
  * Johnston et al. (2020a) Johnston, K. G., Hoare, M. G., Beuther, H., et al. 2020a, ApJ, 896, 35
  * Johnston et al. (2020b) Johnston, K. G., Hoare, M. G., Beuther, H., et al. 2020b, ApJ, 896, 35
  * Johnston et al. (2015) Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19
  * Keto & Wood (2006) Keto, E. & Wood, K. 2006, ApJ, 637, 850
  * Klassen et al. (2014) Klassen, M., Kuiper, R., Pudritz, R. E., et al. 2014, ApJ, 797, 4
  * Klassen et al. (2016) Klassen, M., Pudritz, R. E., Kuiper, R., Peters, T., & Banerjee, R. 2016, ApJ, 823, 28
  * Kobulnicky et al. (2014) Kobulnicky, H. A., Kiminki, D. C., Lundquist, M. J., et al. 2014, ApJS, 213, 34
  * Kolb et al. (2013) Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80
  * Kraus et al. (2017) Kraus, S., Kluska, J., Kreplin, A., et al. 2017, ApJ, 835, L5
  * Krumholz et al. (2007a) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, ApJ, 665, 478
  * Krumholz et al. (2007b) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007b, ApJ, 656, 959
  * Krumholz et al. (2009a) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009a, Science, 323, 754
  * Krumholz et al. (2009b) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009b, Science, 323, 754
  * Kuiper et al. (2010a) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ, 722, 1556
  * Kuiper et al. (2011) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20
  * Kuiper et al. (2010b) Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010b, A&A, 511, A81
  * Kuiper & Yorke (2013a) Kuiper, R. & Yorke, H. W. 2013a, ApJ, 763, 104
  * Kuiper & Yorke (2013b) Kuiper, R. & Yorke, H. W. 2013b, ApJ, 772, 61
  * Laor & Draine (1993) Laor, A. & Draine, B. T. 1993, ApJ, 402, 441
  * Machida et al. (2011) Machida, M. N., ichiro Inutsuka, S., & Matsumoto, T. 2011, The Astrophysical Journal, 729, 42
  * Mahy et al. (2013) Mahy, L., Rauw, G., De Becker, M., Eenens, P., & Flores, C. A. 2013, A&A, 550, A27
  * Maud et al. (2019) Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2019, A&A, 627, L6
  * Maud et al. (2018) Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2018, A&A, 620, A31
  * Maud et al. (2017) Maud, L. T., Hoare, M. G., Galván-Madrid, R., et al. 2017, MNRAS, 467, L120
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Meyer et al. (2019a) Meyer, D. M. A., Haemmerlé, L., & Vorobyov, E. I. 2019a, MNRAS, 484, 2482
  * Meyer et al. (2019b) Meyer, D. M. A., Kreplin, A., Kraus, S., et al. 2019b, MNRAS, 487, 4473
  * Meyer et al. (2018) Meyer, D. M.-A., Kuiper, R., Kley, W., Johnston, K. G., & Vorobyov, E. 2018, MNRAS, 473, 3615
  * Meyer et al. (2021) Meyer, D. M. A., Vorobyov, E. I., Elbakyan, V. G., et al. 2021, MNRAS, 500, 4448
  * Meyer et al. (2022) Meyer, D. M. A., Vorobyov, E. I., Elbakyan, V. G., et al. 2022, MNRAS, 517, 4795
  * Meyer et al. (2019c) Meyer, D. M.-A., Vorobyov, E. I., Elbakyan, V. G., et al. 2019c, MNRAS, 482, 5459
  * Meyer et al. (2017) Meyer, D. M.-A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017, MNRAS, 464, L90
  * Michael & Durisen (2010) Michael, S. & Durisen, R. H. 2010, MNRAS, 406, 279
  * Mignon-Risse et al. (2023a) Mignon-Risse, R., González, M., & Commerçon, B. 2023a, A&A, 673, A134
  * Mignon-Risse et al. (2020) Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42
  * Mignon-Risse et al. (2021) Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2021, A&A, 652, A69
  * Mignon-Risse et al. (2023b) Mignon-Risse, R., Oliva, A., González, M., Kuiper, R., & Commerçon, B. 2023b, A&A, 672, A88
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7
  * Motogi et al. (2017) Motogi, K., Hirota, T., Sorai, K., et al. 2017, ApJ, 849, 23
  * Nayakshin (2010) Nayakshin, S. 2010, MNRAS, 408, L36
  * Nayakshin (2016) Nayakshin, S. 2016, MNRAS, 461, 3194
  * Nayakshin (2017) Nayakshin, S. 2017, PASA, 34, e002
  * Nayakshin & Lodato (2012) Nayakshin, S. & Lodato, G. 2012, MNRAS, 426, 70
  * Oliva & Kuiper (2023a) Oliva, A. & Kuiper, R. 2023a, A&A, 669, A80
  * Oliva & Kuiper (2023b) Oliva, A. & Kuiper, R. 2023b, A&A, 669, A81
  * Oliva & Kuiper (2020) Oliva, G. A. & Kuiper, R. 2020, A&A, 644, A41
  * Papaloizou & Savonije (1991) Papaloizou, J. C. & Savonije, G. J. 1991, MNRAS, 248, 353
  * Persi et al. (1986) Persi, P., Ferrari-Toniolo, M., & Spinoglio, L. 1986, A&A, 157, 29
  * Peters et al. (2010) Peters, T., Banerjee, R., Klessen, R. S., et al. 2010, ApJ, 711, 1017
  * Pickett et al. (2003) Pickett, B. K., Mejía, A. C., Durisen, R. H., et al. 2003, ApJ, 590, 1060
  * Purser et al. (2018) Purser, S. J. D., Lumsden, S. L., Hoare, M. G., & Cunningham, N. 2018, MNRAS, 475, 2
  * Purser et al. (2016) Purser, S. J. D., Lumsden, S. L., Hoare, M. G., et al. 2016, MNRAS, 460, 1039
  * Rafikov (2005) Rafikov, R. R. 2005, ApJ, 621, L69
  * Rafikov (2007) Rafikov, R. R. 2007, ApJ, 662, 642
  * Regály & Vorobyov (2017) Regály, Z. & Vorobyov, E. 2017, A&A, 601, A24
  * Reiter et al. (2017) Reiter, M., Kiminki, M. M., Smith, N., & Bally, J. 2017, MNRAS, 470, 4671
  * Rosdahl & Teyssier (2015) Rosdahl, J. & Teyssier, R. 2015, MNRAS, 449, 4380
  * Rosen et al. (2016) Rosen, A. L., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2016, MNRAS, 463, 2553
  * Samal et al. (2018) Samal, M. R., Chen, W. P., Takami, M., Jose, J., & Froebrich, D. 2018, MNRAS, 477, 4577
  * Seifried et al. (2011) Seifried, D., Banerjee, R., Klessen, R. S., Duffin, D., & Pudritz, R. E. 2011, MNRAS, 417, 1054
  * Seifried et al. (2013) Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320
  * Seifried et al. (2015) Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2015, MNRAS, 446, 2776
  * Seifried et al. (2012) Seifried, D., Pudritz, R. E., Banerjee, R., Duffin, D., & Klessen, R. S. 2012, MNRAS, 422, 347
  * Stecklum et al. (2017) Stecklum, B., Heese, S., Wolf, S., et al. 2017, ArXiv e-prints [arXiv:1712.01451]
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Vaidya et al. (2011) Vaidya, B., Fendt, C., Beuther, H., & Porth, O. 2011, ApJ, 742, 56
  * Vaidya et al. (2018) Vaidya, B., Mignone, A., Bodo, G., Rossi, P., & Massaglia, S. 2018, ApJ, 865, 144
  * Vorobyov (2016) Vorobyov, E. I. 2016, A&A, 590, A115
  * Vorobyov & Basu (2005) Vorobyov, E. I. & Basu, S. 2005, ApJ, 633, L137
  * Vorobyov & Basu (2006) Vorobyov, E. I. & Basu, S. 2006, ApJ, 650, 956
  * Vorobyov & Basu (2010) Vorobyov, E. I. & Basu, S. 2010, ApJ, 719, 1896
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
  * Vorobyov et al. (2013) Vorobyov, E. I., DeSouza, A. L., & Basu, S. 2013, ApJ, 768, 131
  * Vorobyov & Elbakyan (2018) Vorobyov, E. I. & Elbakyan, V. G. 2018, A&A, 618, A7


[datatype=bibtex] [fieldsource=doi,final] [fieldset=url,null]

# A short proof of the strong three dimensional Gaussian product inequality

Ronan HERRY R.H gratefully acknowledges funding from the Centre Henri
Lebesgue. IRMAR, Université de Rennes 1  Dominique MALICET  LAMA, Université
Gustave Eiffel  Guillaume POLY  IRMAR, Université de Rennes 1

###### Abstract

We prove the strong form of the Gaussian product conjecture in dimension
three. Our purely analytical proof simplifies previously known proofs based on
combinatorial methods or computer-assisted methods, and allows us to solve the
case of any triple of even positive integers which remained open so far.

##  1 Introduction and main result

In this note, we prove the following theorem.

######  Theorem 1.1.

Let $(X_{1},X_{2},X_{3})$ be centered real Gaussian vector, and
$p_{1},p_{2},p_{3}\in 2\mathbb{N}$. Then,

(1.1)
$\operatorname{\mathbb{E}}\left[X_{1}^{p_{1}}X_{2}^{p_{2}}X_{3}^{p_{3}}\right]\geq\operatorname{\mathbb{E}}\left[X_{1}^{p_{1}}\right]\operatorname{\mathbb{E}}\left[X_{2}^{p_{2}}\right]\operatorname{\mathbb{E}}\left[X_{3}^{p_{3}}\right],$

with equality if and only if $X_{1},X_{2},X_{3}$ are independent.

Hence, our result completely solves the case $n=3$ of the strong form of the
celebrated _Gaussian product conjecture_. For short, let us introduce the
following notation.

######  Definition 1.2.

We say that $n\in\mathbb{N}$, and $p_{1},\dots,p_{n}\in(0,\infty)$ satisfy the
_Gaussian product inequality_ , and we write
$\mathbf{GPI}_{n}(p_{1},\dots,p_{n})$ provided for all real centered Gaussian
vector $(X_{1},\dots,X_{n})$:

$\operatorname{\mathbb{E}}\left[\prod_{i=1}^{n}\lvert
X_{i}\rvert^{p_{i}}\right]\geq\prod_{i=1}^{n}\operatorname{\mathbb{E}}\left[\lvert
X_{i}\rvert^{p_{i}}\right],$

with equality if and only if $X_{1},\dots,X_{n}$ are independent.

######  Conjecture 1.3.

For all $n\in\mathbb{N}$, and all $p_{1},\dots,p_{n}\in 2\mathbb{N}$,
$\mathbf{GPI}_{n}(p_{1},\dots,p_{n})$ holds.

Despite having received considerable attention, the general case of the
conjecture had, until now, remained wide open. The previous state of the art
regarding the Gaussian product conjecture was the following.

######  Theorem 1.4.

  1. The following cases of Conjecture 1.3 are known.

  2. (a)

For all $p_{1},p_{2}\in 2\mathbb{N}$, $\mathbf{GPI}_{2}(p_{1},p_{2})$.

  3. (b)

([Fre08]) For all $n\in\mathbb{N}$, $\mathbf{GPI}_{n}(2,2,\dots,2)$.

  4. (c)

([LHS20]) For all $p\in 2\mathbb{N}$, $\mathbf{GPI}_{3}(p,p,p)$.

  5. (d)

([RS22a]) For all $p\in 2\mathbb{N}$, $\mathbf{GPI}_{3}(p,6,4)$ and
$\mathbf{GPI}_{2}(p,2,2,2)$.

  6. (e)

([RS22]) For all $p$ and $q\in 2\mathbb{N}$, $\mathbf{GPI}_{3}(2,p,q)$.

The above results are obtained through sophisticated methods. In particular,
[LHS20] relies on a heavily combinatorial approach in connection with the
theory of Gaussian hypergeometric functions; while [RS22a, RS22] is a
computer-assisted method based on the SOS algorithm which provides an explicit
expansion of a positive multivariate polynomial into a sum of squared
quantities. On the contrary, our approach is purely analytical and combines an
optimization procedure through the use of Lagrange multipliers with Gaussian
analysis. Our contribution not only drastically simplifies the proof of the
known cases in dimension three (Theorem 1.4 c and d), but it also enables us
to fully resolve the three dimensional case, that is to say for every choice
of even integer exponents.

######  Remark 1.5.

Let us also mention that, for all $n\in\mathbb{N}$ and $p_{1},\dots,p_{n}\in
2\mathbb{N}$, [Ari98] establishes the _complex_ counterpart of the conjecture;
while [MNPP16] derives a variant of the inequality involving _Hermite
polynomials_. [Wei14] proves $\mathbf{GPI}_{n}(p_{1},\dots,p_{n})$ for all
$n\in\mathbb{N}$ and $p_{1},\dots,p_{n}\in(-1,0)$. Several authors also derive
weaker form of Conjecture 1.3 by considering only Gaussian vectors with
additional assumptions on the covariance matrix. Among others, let us quote
the two recent contributions: [ERR22] for the case of positive correlations,
and [Oui22] for multinomial covariances.

##  2 Proof of the main result

In the rest of the paper, we fix a probability space
$(\Omega,\mathcal{F},\mathbb{P})$ supporting an independent sequence
$(G_{k})_{k\in\mathbb{N}}$ of centered normalized Gaussian variables on
$\Omega$. In the following, $\Sigma=(\sigma_{i,j})$ is a real symmetric non-
negative matrix of size $n$. To $\Sigma$, we associate a centered Gaussian
vector $\vec{X}=(X_{1},\ldots,X_{n})$ with covariance matrix $\Sigma$ by
setting $\vec{X}=\Sigma^{1/2}\vec{G}$ where $\vec{G}=(G_{1},\ldots,G_{n})$.
Let $p_{1},\dots,p_{n}\in 2\mathbb{N}^{*}$, and
$h(x_{1},\dots,x_{n})=x_{1}^{p_{1}}\dots x_{n}^{p_{n}}$. Our strategy consists
in studying the points where the map $\Phi\colon\Sigma\mathchar
14765\relax\operatorname{\mathbb{E}}[h(X_{1},\dots,X_{n})]$ reaches its
minimum. Our argument allows us to characterize those minimal points for $n=2$
or $3$. Using Wick formula [Jan97, Thm. 1.28], it is readily checked that
$\Phi$ is polynomial in the entries of $\Sigma$. We shall need the following
standard lemma. We recall a proof for the sake of self-containedness.

######  Lemma 2.1.

Let $(X_{1},\cdots,X_{n})$ be a Gaussian vector that is centred with
covariance matrix $\Sigma$ non necessarily invertible. Then it holds:

(2.1)
$\displaystyle\operatorname{\mathbb{E}}\left[X_{i}h(X_{1},\ldots,X_{n})\right]=\sum_{j=1}^{n}\sigma_{i,j}\operatorname{\mathbb{E}}\left[\partial_{j}h(X_{1},\ldots,X_{n})\right],\qquad
i\in\\{1,\dots,n\\};$ (2.2)
$\displaystyle\frac{\partial}{\partial\sigma_{i,j}}\operatorname{\mathbb{E}}\left[h(X_{1},\ldots,X_{n})\right]=\operatorname{\mathbb{E}}\left[\partial_{x_{i}}\partial_{x_{j}}h(X_{1},\ldots,X_{n})\right],\qquad
i\neq j\in\\{1,\dots,n\\}.$

###### Proof of Lemma 2.1.

In view of Wick formula, 2.1 and 2.2 are equalities involving polynomials in
the entries of $\Sigma$. It is thus sufficient to establish them for an
invertible $\Sigma$. In this case, let us write $f_{\Sigma}$ for the density
distribution associated with $\vec{X}$. A direct computation yields

$x_{i}f_{\Sigma}+\sum_{j=1}^{n}\sigma_{i,j}\partial_{x_{j}}f_{\Sigma}=0,\qquad
i\in\\{1,\dots,n\\}.$

2.1 readily follows. In order to prove 2.2, consider the Fourier transform of
$f_{\Sigma}$:

$\widehat{f_{\Sigma}}(x_{1},\cdots,x_{n})=\exp\left\lparen-\frac{1}{2}\sum_{i,j=1}^{n}x_{i}x_{j}\sigma_{i,j}\right\rparen,\qquad(x_{1},\dots,x_{n})\in\mathbb{R}^{n}.$

Differentiating this formula, we obtain for $i\neq j$:

$\widehat{\partial_{\sigma_{i,j}}f_{\Sigma}}=\partial_{\sigma_{i,j}}\widehat{f_{\Sigma}}=-x_{i}x_{j}\widehat{f_{\Sigma}}=\widehat{\partial_{x_{i}}\partial_{x_{j}}f_{\Sigma}}.$

Since the Fourier transform is into, we get that
$\partial_{\sigma_{i,j}}f_{\Sigma}=\partial_{x_{i}}\partial_{x_{j}}f_{\Sigma}$,
from which 2.1 readily follows. ∎

In order to highlight the line of reasoning we use in the proof of Theorem
1.1, and for the sake of completeness, let us first give a short proof of
Theorem 1.4 a.

###### Proof of Theorem 1.4 a.

Fix $p_{1}$ and $p_{2}\in 2\mathbb{N}^{*}$. We want to prove that if
$(X_{1},X_{2})$ is a centered real Gaussian vector then
$\operatorname{\mathbb{E}}[X_{1}^{p_{1}}X_{2}^{p_{2}}]\geq\operatorname{\mathbb{E}}[X_{1}^{p_{1}}]\operatorname{\mathbb{E}}[X_{2}^{p_{2}}]$,
with equality if and only if $X_{1}$ and $X_{2}$ are independent. By
homogeneity it is enough to prove the statement when $X_{1}$ and $X_{2}$ are
normalized, in which case the right term of the inequality depends only on
$p_{1}$ and $p_{2}$. Setting, for $t$ in $[-1,1]$,
$\Phi(t)=\operatorname{\mathbb{E}}[X_{1}^{p_{1}}X_{2}^{p_{2}}]$ where
$(X_{1},X_{2})$ is the Gaussian vector associated to

$\Sigma=\begin{pmatrix}1&t\\\ t&1\end{pmatrix},$

the claim is equivalent to the fact that $\Phi$ reaches its unique minimum at
$t=0$. From 2.2,
$\Phi^{\prime}(t)=p_{1}p_{2}\operatorname{\mathbb{E}}[X_{1}^{p_{1}-1}X_{2}^{p_{2}-1}]$
and
$\Phi^{\prime\prime}(t)=p_{1}(p_{1}-1)p_{2}(p_{2}-1)\operatorname{\mathbb{E}}[X_{1}^{p_{1}-2}X_{2}^{p_{2}-2}]$.
In particular, $\Phi^{\prime\prime}>0$ and $\Phi^{\prime}(0)=0$. Consequently,
$0$ is a critical point of a strictly convex function, and thus it is the
unique global minimizer of $\Phi$, from which the result follows. ∎

Theorem 1.1 follows from the recursive argument below; the corresponding
initialization is given by Theorem 1.4 a.

######  Proposition 2.2.

Let $p_{1},p_{2},p_{3}\in 2\mathbb{N}^{*}$. If
$\mathbf{GPI}_{3}(p_{1}-2,p_{2},p_{3})$, then
$\mathbf{GPI}_{3}(p_{1},p_{2},p_{3})$.

###### Proof.

Let $\mathcal{C}$ be the set of real symmetric positive matrices of size $3$
with $1$ on the diagonal, namely

$\mathcal{C}=\left\\{\Sigma=\begin{pmatrix}1&a&b\\\ a&1&c\\\
c&b&1\end{pmatrix}:|a|,|b|,|c|\leq 1,\,\det(\Sigma)\geq 0\right\\}.$

We identify $\mathcal{C}$ with a compact subset of $\mathbb{R}^{3}$. With this
notation, $\mathbf{GPI}_{3}(p_{1},p_{2},p_{3})$ turns out to be equivalent to
the fact that $\Phi$ attains its unique minimum on $\mathcal{C}$ at $I_{3}$.
Since $\mathcal{C}$ is compact and $\Phi$ continuous, $\Phi$ has a global
minimum at some possibly non-unique

$\Sigma_{0}=\begin{pmatrix}1&a&b\\\ a&1&c\\\
c&b&1\end{pmatrix}\in\mathcal{C}.$

We prove that $\Sigma_{0}=I_{3}$. We split the argument in three cases,
depending on the location of $\Sigma_{0}$ in $\mathcal{C}$.

### Case 1.

We assume that $\Sigma_{0}$ is in the interior on $\mathcal{C}$. This means
that $\det(\Sigma_{0})>0$ and $\lvert a\rvert,\lvert b\rvert,\lvert
c\rvert<1$. In this case, $\Sigma_{0}$ is a critical point of $\Phi$. Write

$U=X_{1}^{p_{1}-1}X_{2}^{p_{2}-1}X_{3}^{p_{3}-1}.$

According to 2.2

(2.3)
$\begin{cases}\partial_{a}\Phi(\Sigma_{0})=p_{1}p_{2}\operatorname{\mathbb{E}}[X_{3}U],\\\
\partial_{b}\Phi(\Sigma_{0})=p_{1}p_{3}\operatorname{\mathbb{E}}[X_{2}U],\\\
\partial_{c}\Phi(\Sigma_{0})=p_{2}p_{3}\operatorname{\mathbb{E}}[X_{1}U].\end{cases}$

Thus,
$\operatorname{\mathbb{E}}[X_{1}U]=\operatorname{\mathbb{E}}[X_{2}U]=\operatorname{\mathbb{E}}[X_{3}U]=0$.
On the other hand, let

$V=X_{1}^{p_{1}-1}X_{2}^{p_{2}}X_{3}^{p_{3}}.$

Thus, by 2.1 and the fact that the derivatives vanish,

$\begin{split}\Phi(\Sigma_{0})&=\operatorname{\mathbb{E}}[X_{1}V]\\\
&=(p_{1}-1)\operatorname{\mathbb{E}}[X_{1}^{p_{1}-2}X_{2}^{p_{2}}X_{3}^{p_{3}}]+p_{2}a\operatorname{\mathbb{E}}[X_{3}U]+p_{3}b\operatorname{\mathbb{E}}[X_{2}U]\\\
&=(p_{1}-1)\operatorname{\mathbb{E}}[X_{1}^{p_{1}-2}X_{2}^{p_{2}}X_{3}^{p_{3}}]\end{split}$

In view of $\mathbf{GPI}_{3}(p_{1}-2,p_{2},p_{3})$, we thus get

$\Phi(\Sigma_{0})\geq(p_{1}-1)\operatorname{\mathbb{E}}[X_{1}^{p_{1}-2}]\operatorname{\mathbb{E}}[X_{2}^{p_{2}}]\operatorname{\mathbb{E}}[X_{3}^{p_{3}}]=\operatorname{\mathbb{E}}[X_{1}^{p_{1}}]\operatorname{\mathbb{E}}[X_{2}^{p_{2}}]\operatorname{\mathbb{E}}[X_{3}^{p_{3}}]=\Phi(I_{3}).$

Since $\Sigma_{0}$ is a minimizer, we actually have that
$\Phi(\Sigma_{0})=\Phi(I_{3})$. In particular, this means that we are in the
equality case of $\mathbf{GPI}_{3}(p_{1}-2,p_{2},p_{3})$. If $p_{1}>2$, this
shows mediately that $\Sigma_{0}=I_{3}$. Similarly, if $p_{2}>2$ or $p_{3}>2$,
we conclude in the same way. If $p_{1}=p_{2}=p_{3}=2$, using Theorem 1.4 a, we
deduce that the components of $(X_{1},X_{2},X_{3})$ are pairwise independent.
Since the vector is Gaussian, the conclusion follows.

### Case 2.

We assume that $|a|,|b|,|c|<1$ and $\det(\Sigma_{0})=0$.

$\Sigma_{0}$ is a priori not a critical point of $\Phi$. Since $\Sigma_{0}$ is
a global minimizer on $\mathcal{C}$ it is also a minimizer of $\Phi$ on the
surface

$S=\left\\{\Sigma=\begin{pmatrix}1&a&b\\\ a&1&c\\\
c&b&1\end{pmatrix}:\det(\Sigma)=0\right\\}.$

By the Lagrange multiplier theorem, we conclude that
$\vec{\nabla}\Phi(\Sigma_{0})$ and $\vec{\nabla}(\det)(\Sigma_{0})$ are
colinear (where $\vec{\nabla}=(\frac{\partial}{\partial
a},\frac{\partial}{\partial b},\frac{\partial}{\partial c})$). We have already
computed $\vec{\nabla}\Phi(\Sigma_{0})$ in 2.3, and we have:

######  Lemma 2.3.

$\vec{\nabla}(\det)(\Sigma_{0})=(\alpha_{1}\alpha_{2},\alpha_{1}\alpha_{3},\alpha_{2}\alpha_{3})$,
where $(\alpha_{1},\alpha_{2},\alpha_{3})$ is some non zero vector of
$\ker(\Sigma_{0})$.

###### Proof.

Write $A=2\mathrm{adj}(\Sigma_{0})$ where $\mathrm{adj}$ stands for the
adjugate matrix Since $\det(\Sigma_{0})=0$, $\mathrm{rank}(\Sigma_{0})\leq 2$,
and since $|a|,|b|,|c|<1$, two columns of $\Sigma_{0}$ cannot be proportional
so $\mathrm{rank}(\Sigma_{0})=2$. This implies that $A$ has rank $1$, thus
$A=\alpha^{T}\alpha$ where
$\alpha=(\alpha_{1},\alpha_{2},\alpha_{3})\in\ker(\Sigma_{0})\setminus\\{0\\}$.
By Jacobi’s formula, we have that
$\vec{\nabla}(\det)(\Sigma_{0})=(A_{1,2},A_{1,3},A_{2,3})$. ∎

We deduce that there exists a real number $k$ such that

(2.4) $\begin{cases}\partial_{a}\Phi(\Sigma_{0})=k\alpha_{1}\alpha_{2},\\\
\partial_{b}\Phi(\Sigma_{0})=k\alpha_{1}\alpha_{3},\\\
\partial_{c}\Phi(\Sigma_{0})=k\alpha_{2}\alpha_{3}.\end{cases}.$

Since $(\alpha_{1},\alpha_{2},\alpha_{3})$ belongs to $\ker(\Sigma_{0})$,
$\alpha_{1}X_{1}+\alpha_{2}X_{2}+\alpha_{3}X_{3}=0$ almost surely. From 2.3,

$p_{1}\alpha_{1}\frac{\partial\Phi}{\partial
c}(\Sigma_{0})+p_{2}\alpha_{2}\frac{\partial\Phi}{\partial
b}(\Sigma_{0})+p_{3}\alpha_{3}\frac{\partial\Phi}{\partial
a}(\Sigma_{0})=p_{1}p_{2}p_{3}\operatorname{\mathbb{E}}\left[U\left\lparen\alpha_{1}X_{1}+\alpha_{2}X_{2}+\alpha_{3}X_{3}\right\rparen\right]=0.$

Thus by reporting in 2.4,

$0=(p_{1}+p_{2}+p_{3})k\alpha_{1}\alpha_{2}\alpha_{3}.$

If $k=0$, then 2.4 gives that $\Sigma_{0}$ is a critical point of $\Phi$, and
as in Case 1 we obtain that $\Sigma_{0}=I_{3}$, which contradicts
$\det(\Sigma_{0})=0$. If one of the $\alpha_{i}$ is zero, say $\alpha_{1}$,
then $\alpha_{2}X_{2}+\alpha_{3}X_{3}=0$, so $X_{3}$ and $X_{2}$ are
proportional, and since they are normalized, $X_{3}=\pm X_{2}$, which
contradicts $|c|<1$. Hence, Case 2 cannot happen.

### Case 3.

We assume that $\\{|a|,|b|,|c|\\}\cap\\{1\\}\neq\emptyset$. Say for example
that $|c|=1$. That implies that $X_{3}=\pm X_{2}$ and so by the two
dimensional case Theorem 1.4 a,

$\Phi(\Sigma_{0})=\operatorname{\mathbb{E}}[{X_{1}}^{p_{1}}{X_{2}}^{p_{2}+p_{3}}]\geq\operatorname{\mathbb{E}}[{X_{1}}^{p_{1}}]\operatorname{\mathbb{E}}[{X_{2}}^{p_{2}+p_{3}}]>\operatorname{\mathbb{E}}[{X_{1}}^{p_{1}}]\operatorname{\mathbb{E}}[{X_{2}}^{p_{2}}]\operatorname{\mathbb{E}}[{X_{2}}^{p_{3}}]=\Phi(I_{3}).$

In particular $\Sigma_{0}$ is not a minimizer. It is a contradiction, and Case
3 cannot either happen.

### Conclusion.

We obtain that the only minimizer of $\Phi$ is $I_{3}$ which concludes the
proof. ∎

## References

  * [Ari98] J. Arias-de-Reyna  “Gaussian variables, polynomials and permanents”  In _Linear Algebra Appl._ 285.1-3, 1998, pp. 107–114  DOI: 10.1016/S0024-3795(98)10125-8
  * [ERR22] Dominic Edelmann, Donald Richards and Thomas Royen  “Product Inequalities for Multivariate Gaussian, Gamma, and Positively Upper Orthant Dependent Distributions”  arXiv, 2022  DOI: 10.48550/ARXIV.2204.06220
  * [Fre08] Péter E. Frenkel  “Pfaffians, Hafnians and products of real linear functionals”  In _Math. Res. Lett._ 15.2-3, 2008, pp. 351–358  DOI: 10.4310/MRL.2008.v15.n2.a12
  * [Jan97] Svante Janson  “Gaussian Hilbert Spaces”, Cambridge Tracts in Mathematics  Cambridge University Press, 1997  DOI: 10.1017/CBO9780511526169
  * [LHS20] Guolie Lan, Ze-Chun Hu and Wei Sun  “The three-dimensional Gaussian product inequality” Id/No 123858  In _J. Math. Anal. Appl._ 485.2, 2020, pp. 19  DOI: 10.1016/j.jmaa.2020.123858
  * [MNPP16] Dominique Malicet, Ivan Nourdin, Giovanni Peccati and Guillaume Poly  “Squared chaotic random variables: new moment inequalities with applications”  In _J. Funct. Anal._ 270.2, 2016, pp. 649–670  DOI: 10.1016/j.jfa.2015.10.013
  * [Oui22] Frédéric Ouimet  “The Gaussian product inequality conjecture for multinomial covariances”  arXiv, 2022  DOI: 10.48550/ARXIV.2209.01505
  * [RS22] Oliver Russell and Wei Sun  “Moment ratio inequality of bivariate Gaussian distribution and three-dimensional Gaussian product inequality”  arXiv, 2022  DOI: 10.48550/ARXIV.2208.13957
  * [RS22a] Oliver Russell and Wei Sun  “Using Sums-of-Squares to Prove Gaussian Product Inequalities”  arXiv, 2022  DOI: 10.48550/ARXIV.2205.02127
  * [Wei14] Ang Wei  “Representations of the absolute value function and applications in Gaussian estimates”  In _J. Theor. Probab._ 27.4, 2014, pp. 1059–1070  DOI: 10.1007/s10959-013-0486-z