# Non-LTE Monte Carlo Radiative Transfer. III. The thermal properties of
Tilted and Warped Be Star Discs

M.W. Suffak,1 C.E. Jones,1 A.C. Carciofi,2 T.H. de Amorim2
1Department of Physics and Astronomy, Western University, London, ON N6A 3K7,
Canada
2Instituto de Astronomia, Geofísica e Ciéncias Atmosféricas, Universidade de
São Paulo, Brazil
E-mail<EMAIL_ADDRESS>
(Accepted 2023 September 07. Received 2023 August 28; in original form 2023
June 19 )

###### Abstract

We use the three-dimensional Monte Carlo radiative transfer code hdust to
model Be stars where the disc is tilted from the equatorial plane of the star.
We compute 128 models across 4 spectral types, B0, B2, B5 and B8, tilting the
disc by $None$, $None$, $None$, and $None$, while varying disc density
according to spectral type. We also compute every model for an average and
high stellar rotation rate. We first discuss non-tilted disc temperatures and
show its non-linear dependence on stellar and disc parameters. We find that
tilting the disc minimally affects the density-weighted average disc
temperature, but tilting does create a temperature asymmetry in disc cross
sections, which is more pronounced for a faster rotation rate. We also
investigate the effect tilting has on $V$-band magnitude, polarization, and
the H$\rm\alpha$ line. Tilting the disc does affect these observables, but the
changes are entirely dependent on the position of the observer relative to the
direction of tilt. We find the observables that distinguish tilting from a
change in density or geometry are the H$\rm\alpha$ line shape, where it can
transition between single-peaked and double-peaked, and the polarization
position angle, whose value is dependent on the projected major elongation
axis of the disc on the sky. We also present one early and one late-type model
with warped discs. We find their temperature structure varies a small amount
from the uniformly tilted models, and the different observables correspond to
different tilt angles, consistent with their expected volume of origin within
the disc.

###### keywords:

binaries: general – circumstellar matter – radiative transfer – stars:
emission-line, Be

††pubyear: 2023††pagerange: Non-LTE Monte Carlo Radiative Transfer. III. The
thermal properties of Tilted and Warped Be Star Discs–B

##  1 Introduction

Classical Be stars are defined as non-supergiant B-type stars that have, or
have had, Balmer lines in emission (Collins, 1987). These emission lines are
known to form in a gaseous circumstellar disc that has developed around the
equator of the star. The exact process which leads to the formation of these
discs is uncertain, but coupling rapid rotation with non-radial pulsations
(Baade et al., 2016) is thought to be the most-likely mechanism for stellar
mass-loss. In addition to Balmer line emission, Be star discs are also
characterised by excess continuum emission, particularly at infrared (IR) and
radio wavelengths, and by linear polarization (for recent examples, see
Ghoreyshi et al., 2021; Marr et al., 2022). The most recent comprehensive
review of classical Be stars is given by Rivinius et al. (2013).

Observables seen from Be star discs are not only highly dependent on the
density structure of the disc, but also on the disc temperature, as the
temperature ultimately is what sets the state of the gas through its level
populations and ionization state (Carciofi & Bjorkman, 2008). Until the late
twentieth century, the temperature of Be star discs was assumed to be
constant, or simply fall off with a radial power-law (Waters, 1986). The first
attempt to self-consistently determine the disc temperature was performed by
Millar & Marlborough (1998), who determined the temperature by equating the
rates of energy gain and loss at each point in the disc. They applied this
technique to various case studies of both early and late-type Be stars in
subsequent publications (Millar & Marlborough, 1999b, c; Millar et al., 2000)
and found temperature differences of thousands of Kelvin between the midplane
and upper edge of the disc. Jones et al. (2004) added to this method by
accounting for metals with the inclusion of iron. They found that for the
early type star, $\gamma$ Cas, the inclusion of metals lead to an overall
cooling of the disc, and a slight heating at the inner most disc, within 3
stellar radii. However for the late type star 1 Del, the most heating occurred
on the outer edges of the disc, which was illuminated by light from the poles,
and the greatest cooling happened in the middle portion of the disc, not near
the dense equatorial plane.

Carciofi & Bjorkman (2006) investigated the temperature structure of early-
type Be star discs with their 3-dimensional (3D) non-local thermodynamic
equilibrium (non-LTE) Monte Carlo code hdust. In their models, they found the
temperature at the dense midplane of the disc initially drops within 3-5
stellar radii before rising back to the optically thin equilibrium
temperature, while the thin upper layers of the disc were approximately
isothermal, consistent with Millar & Marlborough (1998). They also found the
disc to be almost completely ionized, except for a small portion in the
midplane near the minimum temperature. Carciofi & Bjorkman (2008) further
investigated this non-isothermal structure by presenting a self-consistent
solution for the viscous decretion disc (VDD) scenario. They determined that
the varying temperature affects the density structure in two ways: 1) the
radial temperature gradient changes the radial fall-off of the density, and 2)
the reduction in temperature within the midplane results in the collapse of
the disc onto itself, thereby causing a decrease in its scale height. They
conclude that a non-isothermal disc density model must be used for detailed
modelling of Be star disc observables. However, many successful modelling
efforts of Be stars using hdust have utilized the simpler isothermal density
formula while solving for non-isothermal disc temperatures (Silaj et al.,
2016; Ghoreyshi et al., 2018; Suffak et al., 2020; Marr et al., 2021).

Over the past decade, the possibility of Be star discs warping, tilting, and
precessing, has gained a lot of attention (see Martin et al., 2011; Brown et
al., 2019; Suffak et al., 2022, for example). There have been a number of
studies using hydrodynamical simulations to predict the nature of warping,
tilting, and oscillations of Be star discs in situations where a binary
companion’s orbit is misaligned to the initial plane of the disc (Martin et
al., 2014b; Cyr et al., 2017; Suffak et al., 2022), many of which focus on
Be/X-ray binary system parameters (Martin et al., 2014a; Brown et al., 2019,
for example). The simulations of Suffak et al. (2022) showed that, under the
influence of a misaligned binary companion, a Be star disc can undergo
episodes of disc tearing, as well as develop eccentric gaps near the primary
star during disc dissipation, in addition to tilting, warping, and precessing.
The phenomena of disc precession and disc tearing are the best current
explanation for the behaviour of the observables in the Be star Pleione (Marr
et al., 2022; Martin & Lepp, 2022).

So far, none of the studies that investigated dynamically simulated disc
tilting, warping, precession, etc., have investigated the effects this would
have on the disc temperature structure, or its observables in a systematic
way. As well, late-type Be stars have been dramatically understudied compared
to their early-type counterparts. In this paper, we first provide results of
static 3D radiative transfer models, showing the temperature structure of non-
tilted Be star discs, ranging in spectral type from B0 to B8 (Section 2). We
then show the same discs, uniformly tilted from the equatorial plane, and
discuss their temperature structure (Section 3), before we present two
scenarios where the continuum, Balmer line, and polarization signatures could
allow a tilted disc to be detected (Section 4). We also briefly discuss how a
warped disc may differ from a flat-tilted disc in Section 5. Our discussion
and conclusions are presented in Section 6.

##  2 Non-tilted Disc Temperatures

We chose a computational grid of four spectral types from B0 to B8, to capture
both early and late-type Be star behaviour. The stellar parameters for each
spectral type were taken from Cox (2000) and Silaj et al. (2010), who
interpolated their parameters from Cox (2000). We model our disc density based
on the widely-used equation

$\rho(r,z)=\rho_{0}\bigg{(}\frac{R_{*}}{r}\bigg{)}^{n}\exp{-\frac{z^{2}}{2H^{2}}},$
(1)

where $\rho_{0}$ is the base density, $R_{*}$ is the equatorial radius of the
star, $H$ is the disc scale height, and $r$ and $z$ are respectively the
radial and vertical coordinates in the disc. Equation 1 is physically
motivated by the viscous decretion disc (VDD) model of Lee et al. (1991), that
in its simplest form predicts $n\,=\,3.5$ for an isothermal and geometrically
thin disc, and has been used in many studies, such as Silaj et al. (2010);
Jones et al. (2008); Suffak et al. (2022); Marr et al. (2021). The scale
height is calculated by

$H(r)=\frac{a}{\Omega}\bigg{(}\frac{r}{R_{*}}\bigg{)}^{1.5},$ (2)

where $a$ is the sound speed, calculated assuming a disc temperature 60% of
the star’s effective temperature (Carciofi & Bjorkman, 2006), and $\Omega$ is
the Keplerian orbital frequency at the equator of the star. We selected two
base density ($\rho_{0}$) values for each spectral type, based on the limits
of base density versus stellar effective temperature shown in figure 8a of
Vieira et al. (2017). Figure 8b of Vieira et al. (2017) also shows there is no
bounds on $n$ with effective temperature, so we choose to use values of $n$ of
2 and 3.5 for every spectral type, as these are approximately the lower and
upper limits of $n$ for the majority of stars studied in Vieira et al. (2017).
Finally, we compute each model for two different stellar rotation rates,
setting the critical fraction, W (defined in equation 6 of Rivinius et al.,
2013, as the ratio of the rotational velocity at the equator to the Keplerian
circular orbital velocity at the equator), to 0.7 or 0.95. Figure 9 of
Rivinius et al. (2013) shows 0.7 to be about the average $W$ for Be stars,
while 0.95 is on the extreme upper end, nearing the critical rotation rate
where the outward centrifugal force at the equator would be equal to the
inward pull of gravity.

The disc size is held constant for each spectral type at 50 equatorial radii
($R_{eq}$). The equatorial radius was scaled to be consistent with the chosen
value of $W$, satisfying the formula (Rímulo et al., 2018)

$W=\sqrt{2\bigg{(}\frac{R_{\text{eq}}}{R_{\text{p}}}-1\bigg{)}},$ (3)

where $R_{\text{p}}$ is the stellar polar radius. Table 1 presents stellar and
disc parameters in our models.

Table 1: Stellar and disc parameters used in our hdust grid of models. Left to right is the spectral type, stellar mass, polar radius, fraction of critical rotation velocity, effective temperature, luminosity, disc base density, disc density slope, and model number. Sp. Type | M ($\rm M_{\odot}$) |  $\rm R_{p}$ ($\rm R_{\odot}$) | $W$ |  $T_{\rm eff}$ (K) |  $L$ ($\rm L_{\odot}$) |  $\rm\rho_{0}$ ($\rm g\,cm^{-3}$) | n | Model #
---|---|---|---|---|---|---|---|---
B0 | 17.5 | 7.4 | 0.7 | 30000 | 39740 | $1\times 10^{-10}$ | 2/3.5 | 1/2
|  |  | 0.7 |  |  | $1\times 10^{-11}$ | 2/3.5 | 3/4
|  |  | 0.95 |  |  | $1\times 10^{-10}$ | 2/3.5 | 5/6
|  |  | 0.95 |  |  | $1\times 10^{-11}$ | 2/3.5 | 7/8
B2 | 9.11 | 5.33 | 0.7 | 21000 | 4950 | $5\times 10^{-11}$ | 2/3.5 | 9/10
|  |  | 0.7 |  |  | $5\times 10^{-12}$ | 2/3.5 | 11/12
|  |  | 0.95 |  |  | $5\times 10^{-11}$ | 2/3.5 | 13/14
|  |  | 0.95 |  |  | $5\times 10^{-12}$ | 2/3.5 | 15/16
B5 | 5.9 | 3.9 | 0.7 | 15000 | 690 | $5\times 10^{-12}$ | 2/3.5 | 17/18
|  |  | 0.7 |  |  | $5\times 10^{-13}$ | 2/3.5 | 19/20
|  |  | 0.95 |  |  | $5\times 10^{-12}$ | 2/3.5 | 21/22
|  |  | 0.95 |  |  | $5\times 10^{-13}$ | 2/3.5 | 23/24
B8 | 3.8 | 3.0 | 0.7 | 12000 | 167 | $1\times 10^{-12}$ | 2/3.5 | 25/26
|  |  | 0.7 |  |  | $1\times 10^{-13}$ | 2/3.5 | 27/28
|  |  | 0.95 |  |  | $1\times 10^{-12}$ | 2/3.5 | 29/30
|  |  | 0.95 |  |  | $1\times 10^{-13}$ | 2/3.5 | 31/32

###  2.1 Azimuthally-Averaged Temperature Slice

Across all of our models, regardless of spectral type, rotation rate, or disc
density parameters, we find the following common traits: (1) the tenuous upper
layers of the disc (i.e., far from the midplane) are fully ionized and
approximately isothermal; (2) the very dense disc midplane contains the
diversity between the models, it can be cooler or hotter than the upper layers
depending on model parameters, and is only partially ionized; (3) between
these two regions exists a transition layer between the fully ionized outer
disc and partly ionized inner disc, where relatively thin, hot sheaths arise.
However, when examining the temperature structure in more detail, the
behaviour seen from one model to another is very non-linear and is coupled to
the disc density structure, spectral type, and stellar rotation rate. We
summarize these particularities in appropriate detail below.

####  2.1.1 Early Spectral Types

Figure 1 shows the azimuthally-averaged cross sections of the non-tilted disc
models for the B0 and B2 models of our grid. We see that the models where
$n=2$, whose discs have a slow density fall-off and thus are much more dense
than $n=3.5$, have a very large cool, partially ionized region surrounding the
midplane of the disc, while the outer regions are much hotter and fully
ionized. The inner cool regions are due to the disc being optically thick,
while in the outer regions, the density drops and the temperature can reach
much higher values. In the models where $n=3.5$, we see that the midplane has
a much smaller cool region and then transitions to hotter temperatures with
increasing radius as the disc becomes optically thin. The upper hot layers are
also larger than in the $n=2$ case due to the densities falling off faster and
more of the disc being optically thin. These inner and outer regions are
separated by a hot thin sheath, which has also been seen in other publications
(Carciofi & Bjorkman, 2006; Sigut & Jones, 2007).

In Figure 2, we have plotted the disc temperature and ionization fraction
(fraction of hydrogen in the disc that is ionized) for a column at a radial
position 30 $R_{*}$ for model 11, which prominently displays these hot
sheaths. We can see that the spike in temperature (i.e., the hot sheath)
occurs right at the boundary between the cooler inner portion where the disc
is partially ionized, and the hotter outer layers where the disc is fully
ionized. This can be explained by the inner cool region being optically thick
and locally trapping the UV radiation, so as the vertical direction offers the
largest escape probability due to lower opacity, the UV radiation travels
vertically and further heats the gas directly above and below the inner cool
region. When we compute the bound-free and bound-bound optical depths, as well
as the hydrogen level populations, they show inverse profiles to the
ionization fraction in Figure 2, being highest around the midplane of the disc
and trending towards zero as height in the disc increases. These same trends
occur for all models that have these hot sheaths in their cross sections,
including both early- and late- spectral types.

The position of these hot sheaths also noticeably changes between models, as
they move closer together as the $n$ value rises, or as the disc base density
decreases. Both of these changes to the density structure make the upper disc
layers more tenuous, which allows the disc to be fully ionized to a greater
vertical depth, making the inner partially ionized disc portion thinner, and
thus the transition regions closer to the midplane. This can easily be seen by
comparing the cross sections of model 1 to model 2 and 3, respectively. There
are also cases where the disc is so tenuous that it is nearly entirely
ionized, even in the midplane, seen for example in model 4. Here there is no
cooler section in the midplane. Instead the midplane is hotter than the upper
disc layers, due to the denser midplane being able to reprocess UV radiation
and increase the role of diffuse radiation in that area.

In models 5 to 8 and 13 to 16, we see that increasing the stellar rotation
rate to $W=0.95$ does not change the qualitative temperature patterns in the
disc. However the temperatures in the upper disc are notably higher than in
the slower rotating case and the hot sheaths and disc midplane can be slightly
warmer than with slower rotation. This indicates that the hotter stellar poles
caused by increased gravity darkening at this high rotation is able to “carve"
into the disc deeper, penetrating the disc midplane with more UV radiation and
raising its temperature.

Figure 1: Disc cross sections of the azimuthally-averaged temperatures of
models 1-16, as noted in the subplot titles. The colour of each cell in the
cross section corresponds to temperature as shown in the colour bar. The axes
are in stellar equatorial radii, which is different for each spectral type,
consistent with the polar radius listed in Table 1. Figure 2: Ionization
fraction (blue, left $y$-axis), and temperature (orange, right $y$-axis)
versus $z$ height in the non-tilted disc of model 11. The radial distance of
the measurements from the central star is about 30 $R_{*}$.

####  2.1.2 Late Spectral Types

The cross sections for our B5 and B8 non-tilted models are shown in Figure 3.
In these later spectral types, we see some different behaviour in the
temperature structure than the early B0 and B2 type stars. The highest
density, lowest $n$ model for the B5 spectral type is similar to the analogues
for B0 and B2 stars. However, at lower densities, for $n\,=\,2$ in both B5 and
B8 type stars, we see the midplane is hotter than the outer disc, which is
opposite the early type stars. This is the same as Millar & Marlborough
(1999b) found in their work for the late-type star 1 Del. They explain that
this temperature inversion is due to collisions populating the upper levels,
and thus photoionization from these upper levels is able to heat the gas,
while the disc remains optically thick to Lyman continuum photons. However, as
these hot midplane sections are radially extended in our discs when the
density has already fallen off exponentially, collisions are not going to be a
major factor, and thus this hotter midplane would be due to the discs ability
to reprocess the UV radiation (as mentioned for model 4 in Section 2.1.1) and
locally heat the denser midplane through diffuse radiation. Conversely, we see
in models where $n=3.5$, that the midplane temperature drops off in the
outermost disc. Due to the much faster drop off of density compared to the
$n=2$ model, there is not enough diffuse radiation contributed from the disc
itself to make up for the lack of UV radiation reaching the outer midplane of
the disc from the late-type star.

With a higher rotation rate, we see the B5 models largely retain the same
structure as their slower-rotating counterparts, however the rapidly rotating
B8 models, particularly models 30, 31, and 32, display dramatically hotter
disc midplanes, as well as hot sheaths which did not appear in the slower
rotating case. We interpret this again as the hotter poles being able to
“carve" farther into the disc and cause greater heating in the midplane. Here
the high stellar rotation gives qualities of both an early-type star from the
hot poles, and a late-type star from the very cool equator causing this large
change in temperature cross section.

Figure 3: Same as Figure 1, but for models 17-32. Note the change in the
maximum temperature of the colour bar scale.

##  3 Tilted-Disc Temperatures

(a)

(b)

Figure 4: Two schematics showing the orientation of the tilted disc with
respect to the $x$, $y$, and $z$ axes. The azimuthal angle ($\phi$) is defined
with $\phi=$$$ along the positive $x$-axis as shown. The disc is tilted about
the $x$-axis, and thus the tilt angle $\alpha$ is defined from the $y$-axis.
The central star is represented by the blue circle at the origin of both
diagrams. The disc and star sizes are not to scale.

To expand our work on our non-tilted grid, we tilted all of our models listed
in Table 1 by $\alpha\,=\,$ $None$, $None$, and $None$ away from the
equatorial plane. Figure 4 shows a schematic of the orientation of our tilted
disc models in cartesian coordinates, where the star lies at the origin. The
disc is tilted about the $x$-axis, so the disc tilt angle is measured from the
$y$-axis. When referring to azimuthal angles ($\phi$), we denote the positive
$x$-axis as having $\phi=$$$, and the positive $y$-axis having $\phi=$$$. Thus
$\phi$ values of $None$ and $None$ correspond to the negative $x$ and
$y$-axes, respectively. The height of the disc midplane can thus be given as

$Z(r_{i},\phi_{j})=-r_{i}\sin\bigg{\\{}\arctan\bigg{[}\sin(\phi_{j})\tan\big{(}\alpha\big{)}\bigg{]}\bigg{\\}},$
(4)

where $r_{i}$ and $\phi_{j}$ are the radial and azimuthal coordinates on the
midplane, and $\alpha$ is the disc tilt angle about the $x$-axis (either
$None$, $None$, or $None$).

To assess any global changes in disc temperature due to disc tilting, we
calculate the mass-averaged disc temperature, $\bar{T}_{M}$, using the formula

$\bar{T}_{M}=\frac{1}{M_{\rm disc}}\sum_{i=0}^{N}T_{i}\rho_{i}V_{i},$ (5)

where $M_{disc}$ is the total mass of the disc, and $T_{i}$, $\rho_{i}$ and
$V_{i}$ are the temperature, density, and volume of the $i$-th cell. The sum
is performed over all $N$ cells in the disc. This formula was adapted from the
equation for density-weighted average temperature of McGill et al. (2013). The
results of these calculations are presented in Figure 5. We find that there is
a clear relationship between $\bar{T}_{M}$ and disc tilt angle: the greater
the tilt angle, the greater $\bar{T}_{M}$. While this trend is clear, the
changes are not large, with most discs varying in temperature by less than
1000 K. The change is especially small in our densest models, where the change
in average temperature cannot be distinguished between being caused by the
tilt or simply the error in the nature of Monte-Carlo simulations. We also
note that most of these density-weighted average temperatures are below
$60\%\,T_{\rm eff}$, which Carciofi & Bjorkman (2006) found to be a good
isothermal approximation of Be star discs. Since the density-weighted average
weights the most dense, and hence most optically thick regions the highest, it
is not surprising that our least dense models come closest to, or sometimes
end up greater than, the $60\%\,T_{\rm eff}$ mark. This is consistent with the
results of our non-tilted discs seen in Figures 1 and 3 where the lower
density discs have a higher temperature. It is worth mentioning that by
definition, our tilted models are not anchored at the equator of the star,
thus the innermost disc may have slightly inflated temperatures due to
directly seeing a hotter part of the star than if it were aligned at the
equator.

(a)

(b)

Figure 5: Density weighted average disc temperatures for our early-type (a)
and late-type (b) models. The model number, listed in Table 1 is on the
x-axis. The tilt angle of each disc model is indicated by the legend. The grey
lines indicate 60% of the star’s $T_{\rm eff}$ for each model.

###  3.1 Detailed Temperature Structure

Figure 6: Temperature structure of the line of nodes cross section (left
column) and the cross section farthest from the equator (right column) of the
tilted discs with parameters of model 1 (see Table 1). The top row is for a
$None$ tilt, middle row for $None$ tilt, and bottom row for a $None$ tilt, as
indicated by the $\alpha$ value in the leftmost plot of each row. The colour
corresponds to the temperature of a given cell, as indicated by the colour bar
on the right.

Although the density-weighted average disc temperature doesn’t change
appreciably with tilt angle, we find that a tilted disc can have significant
changes in the temperature structure of certain areas of the disc. An example
is shown in Figure 6, where we show cross sections of discs with parameters of
model 1, tilted by $None$, $None$, and $None$. In the line of nodes cross
section of the disc ($\phi\,=\,$$$, left column) we see that changes in the
temperature structure from the non-tilted model (top left panel of Figure 1)
are hard to detect. However, in the cross section farthest from the equator
($\phi\,=\,$$$, right column of Figure 6), we see noticeable temperature
difference between the top and bottom of the disc, with the top of the disc
becoming cooler and the bottom of the disc becoming hotter. This is due to the
effect of gravity darkening arising from rapid rotation, resulting in the
stellar equator being cooler than the pole. Thus, the part of the disc that
moves closer to the equator when tilted, ends up being cooler than the non-
tilted solution, and the part that moves closer to the stellar pole ends up
hotter than the non-tilted case. We note that the disc midplane does not
significantly change temperature when tilted because its high density compared
to the rest of the disc makes it insensitive to the change in radiation input
caused by the tilting of the disc.

Figure 7: Same as Figure 6, but for the parameters of model 30. Note the
change in the maximum temperature of the colour bar scale.

In Figure 7, we show the same plots as Figure 6, but for one of our late-type
models. We see that the temperature trends of the $\phi\,=\,$$$ cross section
of the disc are the same, with the top of the disc cooling and the bottom
heating as it is oriented closer to the pole of the star. In the line of nodes
of the disc, however, we see a change in the temperature structure: the hot
bands that were on either side of the midplane in the non-tilted models,
greatly lessen in temperature as the overall disc tilt angle increases. Since
the central star does not change, this change in temperature structure is due
to the diffuse radiation field of the disc affecting the disc temperature
significantly. Plots similar to Figures 6 and 7 for all other models are
presented in Appendix A.

(a)

(b)

Figure 8: Plots of temperature vs. radius of model 4 at the midplane of the
disc in two different azimuthal directions; $\phi\,=\,$$$ (a), and
$\phi\,=\,$$$ (b). The four different lines are for different disc tilt angles
as indicated by the legend.

Overall we see that the average cell temperature of the $\phi\,=\,$$$ cross
sections of the disc can differ by as much as 30% when compared to the non-
tilted models. The large majority of this difference however, is in the
optically thin outer disc. As shown in Figure 8, while the midplane
temperature does differ in the $\phi\,=\,$$$ cross sections of the disc from
the non-tilted model, the difference is quite small. Appendix B contains
similar plots of the midplane for all models. As the tilt angle increases, we
see the temperature of the innermost disc increase as well. This is due to the
midplane being oriented closer to the pole, and therefore directly seeing a
hotter area of the star. We also see that the midplane temperature profile of
the tilted disc still has the same structure as past publications (Millar &
Marlborough, 1998; Carciofi & Bjorkman, 2006), with the temperature reaching a
minimum within the first few stellar radii before increasing to an
approximately isothermal temperature in the outer disc. It is important to
note that for some of our densest models, the midplane temperature does not
increase from its minimum due to the high density of the disc. It is also
noteworthy that this structure becomes less pronounced as one moves towards
the top and bottom of the disc, away from the midplane.

##  4 Tilted Disc Observables

Due to the 3-dimensional nature of our simulations, we specify our observer
position by two spherical coordinates, the polar angle $\theta$ and the
azimuthal angle $\phi$. $\theta$ ranges from $None$ when looking pole-on with
the star, to $None$ when looking at the equatorial plane, while $\phi$ is
defined in the range [$None$,$None$) in the same manner as shown in Figure 4.

Observables of Be stars are highly dependent on the orientation of the disc
with respect to the observer. It is well known that a non-tilted disc seen
edge-on will have dimmer photometric magnitudes, higher polarization levels,
and shell emission lines, compared to a disc viewed pole-on, which will have
brighter photometric magnitudes, close to zero polarization, and single-peaked
emission lines. However, in a tilted disc scenario, looking edge-on with the
disc may not imply one is looking equator-on with the star, and a face-on view
of the disc would not be pole-on with the star. Thus, the changes in
observables that occur due to disc tilting would be expected to be highly
dependent on the orientation of the tilt with respect to the observer.

In their hydrodynamic simulations, Suffak et al. (2022) showed that a Be star
disc can tilt and then precess under the influence of a misaligned binary
companion. We now examine these two scenarios where a tilted disc may be able
to be detected through a change in their observables.

###  4.1 Case I: Viewing a Disc with Varying Tilt Angle

In this first scenario, a misaligned binary companion torques the disc above
and below the equatorial plane of the primary star, meaning the axis of the
disc’s tilt is aligned with the companion’s line of nodes. Here, the
orientation between the disc tilt axis and the stationary observer is
constant, thus the observer merely sees the disc tilt in a constant direction
over time.

Since our observing coordinates are defined with respect to the central star
and tilt axis of the disc, we can simulate a disc tilting over time by
plotting any observable viewed from the same ($\theta$, $\phi$) observer
position over four simulations where a model disc is tilted by $None$, $None$,
$None$, and $None$. An example of this is shown in Figure 9, which shows, for
a constant $\theta$ of $None$ and $\phi$ of $None$, $None$, and $None$,
$V$-band magnitude, H$\alpha$ equivalent width (EW), H$\alpha$ violet-to-red
(V/R) ratio, polarization percentage and position angle (PA) in the $V$-band,
versus disc tilt angle for systems with parameters of model 1. Here we see the
change in an observable as the disc tilts is very dependent on the observer’s
position.

For example, at $\phi\,=\,$$$ or $None$, the observer will be aligned with the
tilt axis of the disc, thus the $V$ magnitude, EW, and percent polarization do
not vary as much as at other $\phi$ angles, however the polarization position
angle will vary greatly with increasing tilt angle. This is in contrast to
$\phi\,=\,$$$, where the disc is tilted to be more face-on with the observer.
As the tilt angle increases, we see $V$ magnitude increase, while EW and
polarization decrease and the position angle stays constant. Not plotted is
$\phi\,=\,$$$, where the disc tilts to be more edge-on with the observer, and
the trends would be generally reversed from the $\phi\,=\,$$$ case. Finally,
for an intermediate azimuthal angle of $None$, all of the observables vary
greatly, with the only clear trend being the polarization position angle. Note
that these results are degenerate with $\phi\,=\,$$$, while the degenerate
pair of $None$ and $None$ will give similar results. The polarization position
angle is thus the one observable which shows constant change as the disc tilts
(except when the observer is exactly perpendicular to the tilt axis). We also
see for this simulation, large V/R ratios of about 10% at certain azimuthal
observing angles when the disc is tilted by $None$. This is an extreme example
however, since most of our models don’t reach a V/R ratio of more than 2 to
5%.

Figure 9: Plots showing (top to bottom) $V$-band magnitude, H$\alpha$
equivalent width, H$\alpha$ V/R ratio, polarization percentage, and
polarization position angle (PA) in the $V$-band, versus disc tilt angle, for
systems with parameters of model 1 (Table 1). All points are for a $\theta$
observing angle of $None$, and the different coloured lines indicate different
$\phi$ observing angles as indicated by the legend. Some $\phi$ directions may
be degenerate and thus not every line will show on every plot. Figure 10: Same
format as Figure 9, but for parameters of model 29.

Figure 10 shows the same plot as Figure 9, but for a late-type B8 star, model
29. We see another example of the variability of the change of observables
depending on observer position, and that the magnitude of these changes for
late-type stars are much less than those of the early-type example shown
previously, particularly the $V$ magnitude and polarization, which is expected
due to the lower disc density.

To briefly investigate disc obscuration effects, we took two models, 10 and
25, at a tilt angle of $None$ and computed them with a disc radius of 20
$R_{*}$ instead of 50. These models presented no difference in observables
aside from a slightly lower H$\rm\alpha$ EW. This is expected, as the
H$\rm\alpha$ emission comes from a large portion of the disc, while visible
continuum emission and polarization come from the inner few stellar radii.
This shows that the material in the outer disc is largely optically thin, and
thus the size of a tilted disc would not affect observables other than
H$\rm\alpha$ until the disc radius is reduced to a few stellar radii.

####  4.1.1 H$\rm\alpha$ Line Shapes

The shape of the H$\rm\alpha$ line in Be stars is largely due to Doppler shift
caused by the relative velocity between the disc material and the observer. In
a non-tilted disc, the H$\rm\alpha$ line in Be stars is seen as single-peaked
when looked pole-on with the star and face-on with the disc ($\theta\,=\,$$$),
and is a double-peaked line when observed at other inclination angles due to
disc material moving both toward and away from the observer. As the
H$\rm\alpha$ emission line is a defining characteristic of Be stars, we also
looked at the shapes of the emission lines in our simulations to see if the
changing shape might show indications of disc tilting.

Overall we find there are three different patterns in which the tilt of the
disc can affect the H$\rm\alpha$ line. The first is shown in Figure 11, where
$\phi\,=\,$$$ and the observer is aligned with the tilt axis of the disc, so
the disc is seen as rotating either clockwise or counter-clockwise from the
observers perspective. In this case, we see the line strength stay
approximately the same, but the V/R ratio increases with increasing tilt
angle. This effect is particularly noticeable in the early-type stars at high
($\theta\,\geq\,$$$) inclinations, where the projected area of the disc does
not change very much with increasing tilt angle. The lines plotted here are
one of the most extreme V/R ratios that we have obtained from our models.

The second pattern is shown in Figure 12, where $\phi\,=\,$$$ and the disc
tilts to be more face-on or edge-on with the observer. Here we see that the
line starts out as double-peaked for $None$ tilt, and transitions to a single-
peaked line when the tilt is $None$. This behaviour occurs for all spectral
types and densities. The reverse of this process is also seen in our models,
with the lines shifting from single-peaked to double-peaked with increasing
tilt angle for pole-on observing angles. It is worth noting that the
equivalent width of the line is approximately constant in these scenarios as
well.

Finally when $\phi\,=\,$$$, the motion of the disc relative to the observer is
a combination of a rotation clockwise or counterclockwise in the plane of the
sky, and a tilt to be more edge-on or face-on. Figure 13 shows this case,
where the line shape and peak separation stay the same, but the line decreases
in strength as tilt angle increases for the B0 and B2 type stars, and slightly
increases in strength for the late type stars. The reverse is also seen for
some observing angles. These changes in the normalized lines are largely due
to a change in continuum level rather than a change in the emission itself.

Figure 11: Simulated H$\rm\alpha$ lines of models 1 (top left), 9 (top right),
17 (bottom left) and 25 (bottom right), for four different disc tilt angles as
indicated by the legend. The model spectra are seen from an observer at
position $\phi=$$$, $\theta=$$$. Figure 12: Simulated H$\rm\alpha$ lines of
models 4 (top left), 12 (top right), 20 (bottom left) and 28 (bottom right),
for four different disc tilt angles as indicated by the legend. The model
spectra are seen from an observer at position $\phi=$$$, $\theta=$$$. Figure
13: Simulated H$\rm\alpha$ lines of models 6 (top left), 14 (top right), 22
(bottom left) and 30 (bottom right), for four different disc tilt angles as
indicated by the legend. The model spectra are seen from an observer at
position $\phi=$$$, $\theta=$$$.

###  4.2 Case II: Watching a Tilted Disc Precess

The second scenario where a tilted disc may be able to be detected is where an
already tilted disc precesses under the influence of a misaligned binary
companion. This occurs in many simulations of Suffak et al. (2022),
particularly after mass-loss from the disc is turned off, the disc is no
longer anchored to the equator of the star, and the line of nodes of the disc
is free to rotate about the primary star.

By holding the polar viewing angle ($\theta$) constant, and moving around the
star and disc in $\phi$, we can see what observational signature a disc may
present if it were precessing about the pole of the star. This is shown in
Figure 14, where we’ve plotted the same quantities as Figure 9 versus $\phi$,
for model 1 with a disc tilt of $None$. From this Figure we can see that, as
the observational viewpoint moves around the star/disc system at a constant
$\theta$, the observables oscillate quite significantly as the disc moves from
being edge-on with the observer at $None$ and $None$, to being more face-on at
$None$ and $None$. Moving the observer like this is exactly the same as if the
disc were rigidly precessing about the pole of the star and the observer was
stationary, but this is accomplished here without the need for computationally
expensive hydrodynamical simulations.

Unlike the case of tilting a disc, precession shows signatures in all
observables. The percent polarization and $V$ magnitude will oscillate at half
the precession period as more or less of the inner disc becomes visible as the
disc precesses about the star. The H$\rm\alpha$ EW also oscillates at half the
precession period and can increase with $V$ magnitude as more of the disc is
visible, or decrease as $V$ magnitude increases due to the changing continuum
level. The position angle will oscillate about zero at the same period as the
precession. These period/half-period trends can be seen easily with the dashed
line in Figure 14 shifted $None$. We note here that the half-period trends of
$V$ magnitude, H$\rm\alpha$ EW, and percent polarization are not perfectly
symmetric because the $\theta$ of $None$ means the observer will see slightly
more or less of the central star from opposite sides of the disc. If $\theta$
were $None$, the half-period trends would ideally have perfect symmetry. The
V/R trend is more complex, as the trend from $\phi\,=\,$$$ to $None$ is the
reverse of the trend from $None$ to $None$. This is due to the relative
velocities “reversing" as the observer moves to the opposite side of the disc,
thus the trend in this case is antisymmetric, with $None$ and $None$ being the
nodes of the oscillation.

Figure 14: Top to bottom, the $V$-band magnitude, H$\alpha$ equivalent width,
H$\alpha$ V/R ratio, polarization percentage in the $V$-band, and position
angle, versus azimuthal viewing angle $\phi$, for model 1 with the disc tilted
$None$. The system is viewed at a $\theta$ of $None$. The dashed line is
shifted by $None$ to facilitate comparison between periods.

##  5 Warped vs. Tilted Discs

We recognize our flat tilted models are limited, particularly at the star-disc
boundary, where we have the inner disc tilting the same amount as the outer
disc. In reality, it is more likely that the Be star disc would be anchored to
the equator of the star, and the rest of the disc would be warped away from
the equator by some degree.

To test the difference this may cause, we chose to apply a warp to models 10
and 25, instead of a flat tilt. To warp the computational grid, we fix the
first radial bin at the equator, and then linearly increment the degree of
tilt with each subsequent radial bin, up to a maximum tilt of $None$ at the
furthest radius. The height of the disc midplane is then given by

$Z(r_{i},\phi_{j})=-r_{i}\sin\bigg{\\{}\arctan\bigg{[}\sin(\phi_{j})\tan\big{(}\alpha\frac{i}{49}\big{)}\bigg{]}\bigg{\\}},$
(6)

where the only difference from Equation 4 is that now we have altered the disc
tilt angle $\alpha$, such that it gets bigger with increasing radius ($i$
denotes the radial cell index, which goes from 0-49 in our models), so the
disc becomes warped instead of a flat tilt. These models also have a disc
radius of 50 $R_{*}$, and thus the most highly inclined outer parts of the
disc do not contribute much to the observables in question here as they are
known to originate in the inner disc, which is only moderately tilted.

(a) Model 10, full view

(b) Model 10, zoomed in to 10 $R_{*}$.

(c) Model 25, full view

(d) Model 25, zoomed in to 10 $R_{*}$.

Figure 15: Each panel is a similar to Figure 6; (a) and (b) are for model 10,
while (c) and (d) are for model 25. The top row in each subplot is for the
model with a flat $None$ tilt, while the bottom row is for the model warped to
a maximum of $None$.

Figure 15 shows temperature cross sections of the warped disc compared to its
flat tilted counterpart. We see that in the inner warped disc, the upper and
lower edges of the disc are slightly cooler than the flat tilted model. This
is due to the inner disc still being anchored at the equator in the warped
case, and thus the inner disc is seeing less radiation than when it is tilted
at $None$ with the rest of the disc. This effect is seen even in the non-
warped slice of the disc at $\phi\,=\,$$$, due to less diffuse radiation
within the disc being able to heat this cross section. In the outer warped
disc, the upper and lower edges are respectively cooler and warmer than the
flat-tilted case due to the upper edge being shielded by the inner disc at the
equator, and the lower edge being freely exposed to radiation from higher
stellar latitudes. Figure 15 also shows how, not surprisingly, the temperature
structure warps with the warped density of the disc, and that the interior
structure is the same as the tilted model aside from this warp.

With respect to observables, comparison of the warped disc to the different
tilted simulations highlights how different areas of the disc are responsible
for different observables. To compare these observables between the tilted
models and warped models, we held the $\theta$ observing angle constant, and
calculated a chi-squared value over all $\phi$ observing angles between each
tilted model and the warped model. For the early-type warped model, we find
that the $V$ magnitude, percent polarization, and polarization PA of the
warped model are best matched by the non-tilted and $None$ tilted models, and
occasionally the $None$ tilted model. This is expected, as these observables
originate in the inner disc, which is the least tilted part of the warped
disc. On the other hand, H$\rm\alpha$ EW is best matched by the $None$ tilt
for near pole-on $\theta$ angles, and by $None$ and $None$ tilt angles for
$\theta$ values greater than $None$. This is due to a large increase in the
continuum emission at certain $\phi$ angles for the $None$ flat tilted model,
which causes the H$\rm\alpha$ EW to drop significantly. The V/R ratio of the
H$\rm\alpha$ line also follows the same trends as the EW. The same trends are
seen for the late-type (model 25) warped disc, except that the H$\rm\alpha$
line best matches the flat $None$ tilted disc for all viewing inclinations due
to the optically thin disc.

Figure 16: Same as Figure 14, but for the warped disc of model 10. The system
is viewed at $\theta\,=\,$$$.

For comparison purposes, in Figure 16 we show the observables of a precessing
warped disc, similar to what was shown in Figure 14 for the flat tilted model.
The stellar and disc parameters of both Figures are different, so they cannot
be directly compared, however we do see that a warped disc, if it were
precessing, produces the same period and half-period trends as the flat tilted
model, with some asymmetry in the photometry, polarization, and H$\rm\alpha$
line.

##  6 Discussion and Conclusions

In this paper, we have shown how tilting a Be star disc out of the equatorial
plane of the primary star can affect the disc’s temperature structure, as well
as it’s observables. We modelled B0, B2, B5, and B8 type stars, with different
densities, two different rotation rates, and disc tilt angles of $None$,
$None$, $None$ and $None$.

We find that the temperature structure between non-tilted early and late type
stars can differ greatly, and the behaviour from model to model is highly non-
linear. The exact temperature structure is dependent on the disc density
configuration, the spectral type, and the stellar rotation rate, which means
depending on the model parameters, we see particular trends in the temperature
behaviour.

In our non-tilted models we see that all discs have an inner cool region, the
extent of which dramatically depends on the density exponent $n$. This can be
of significance, as low excitation lines, particularly Fe ii, are known to
originate in these cool inner disc volumes (Carciofi & Bjorkman, 2006). Fe ii
emission lines have been well documented in Be stars (Hanuschik et al., 1996),
and their line-cooling effects have also been explored (Jones et al., 2004).
Since Fe ii emission lines originate in these cool areas, their shape could be
used as a tracer of the radial extent of these regions, assuming the width of
the line is largely due to Doppler broadening. In this sense, if the Fe ii
line had large peak separation and a sharp drop-off in the wings, the center
cool region would be relatively small, however a large cool region would mean
a large formation loci for Fe ii and its line shape could be similar to Balmer
emission lines, albeit with lower peak intensity. Thus, Fe ii lines may be a
valuable constraint on the value of $n$ in Be star disc models, and shows the
great importance of having a non-isothermal disc model (which was attempted
but not conclusively shown by Klement et al., 2015).

We find presence of hot bands above and below the midplane in nearly all disc
density configurations, consistent with findings in other studies (Millar &
Marlborough, 1998; Carciofi & Bjorkman, 2006; Sigut & Jones, 2007; McGill et
al., 2011). We offer the first concrete explanation of these sheaths, showing
in Figure 2 that the sheaths occur right at the boundary between where the
disc is partially and fully ionized. This strongly indicates that these
sheaths are the result of UV radiation that has been trapped in the optically
thick, cold inner disc, escaping vertically through the disc and adding excess
heat to the optically thin outer disc right at the boundary of this partially
ionized region. We also predict that if the disc is dense enough, diffuse
radiation near the midplane of the disc can play a large role in heating the
disc midplane, sometimes causing the midplane to be warmer than the upper disc
layers despite not being fully ionized, particularly seen in our models with
late-type stars.

We also investigated the difference between the star having an average
rotation rate of 70% the critical velocity versus a high rotation rate, 95% of
the critical velocity. In our B0, B2, and B5 models, this increase in rotation
had marginal effects, only heating the outer disc and midplane slightly, but
keeping the overall temperature structure the same. This is different when
looking at our B8 models, where an increase in rotation caused a large
increase in midplane temperature, as well as the appearance of prominent hot
sheaths. In this case, the combination of the higher rotation giving hotter
poles, along with the lower densities used for our B8 models, allow the
radiation from hot poles to “carve" farther into the disc, causing substantial
heating in and around the midplane.

Overall the temperature structure of our non-tilted models is remarkably
similar to the works of Millar & Marlborough (1998, 1999a, 1999b, 1999c);
Millar et al. (2000), who did detailed work on both the early-type star
$\gamma$ Cas, and the late-type star 1 Del, using an escape probability method
and by balancing the energy contributions to calculate the disc temperature
structures of these stars. This similarity comes despite their code using a
different density prescription, and only having five hydrogen energy levels
included. hdust, used here, uses a Monte Carlo technique to solve the
radiative transfer, has 12 non-LTE and 25 LTE hydrogen energy levels, and also
accounts for line radiation and bound-bound transitions. We argue that the
strong agreement of the temperature distributions between previous work,
including Carciofi & Bjorkman (2006) for B3 spectral types also using hdust,
and the new work presented here, infers that the temperature and ionization
levels are primarily controlled by photoionization — recombination
equilibrium. This agreement also provides strong evidence of the broad
applicability of our work to gaseous discs.

In tilting the disc, we see modest large-scale changes to the disc’s average
temperature, with it increasing slightly as the disc tilt angle increases. On
a smaller scale, we find that with increasing tilt angle, the part of the disc
that moves towards the equator becomes cooler, while the portion of the disc
moving towards the stellar pole becomes hotter. These changes are explained by
gravity darkening of the rapidly rotating star, making the stellar poles
hotter than the stellar equator. This anti-symmetric change of disc
temperatures is why the overall average disc temperature does not change
appreciably when the disc is tilted. The temperature in the midplane of the
disc is largely unchanged by the disc tilting due to it’s higher density than
the rest of the disc, although for already optically thin discs, the midplane
can vary in temperature as well. This behaviour is seen across all spectral
types.

Examining the observables of our tilted disc simulations, we offer two
scenarios where a disc tilt may be detected.

The first case is where the disc is actively observed to be tilting. In this
scenario, the change in observables is entirely dependent on the direction of
tilt relative to the observer. A disc may appear completely differently with a
$None$ or $None$ change of relative orientation as the disc either moves to be
more face-on or more edge-on with the observer as it tilts. This variability
would make it difficult to interpret whether changes in observables of a Be
star would be due to a disc tilting or simple changes in disc density or size.
The strongest evidence of disc tilting would appear in the polarization PA
where, if one were looking along the axis of the disc’s tilt, the position
angle should exactly match the tilt angle of the disc, however the position
angle will change some amount as long as the observers line of sight is not
exactly perpendicular to the tilt axis. No other change in geometry would
cause a change of tens of degrees in the position angle, making it a key
measurement to look at for proof of disc tilting. The shape of the
H$\rm\alpha$ line would also be a clear indication of disc tilting as the
shape changes from single-peaked to double-peaked and vice-versa. This change
could not be brought about by a simple change of density structure or a
larger/smaller disc, and could only occur with a major change of disc geometry
such as tilting of the H$\rm\alpha$ emitting region. This would also be seen
in other emission lines as well, not just H$\rm\alpha$. The advantage of these
two observables being the leading indicators of tilting is that one of them
should appear no matter the orientation of the disc, given a large enough disc
tilt. If the disc is tilting to be more face-on or edge-on with the observer,
the H$\rm\alpha$ line shape would change in shape while the position angle
would not. On the other hand, if the observer was looking more aligned with
the tilt axis of the disc, the polarization position angle would change, while
the H$\rm\alpha$ line shape would be approximately constant. Thus, both the
emission line shape and polarization position angle are key signatures of disc
tilting. Another difference that could set a tilted disc apart from a non-
tilted disc is the V/R ratio in the H$\alpha$ emission lines, however this
ratio is not particularly strong in our models apart from a few cases of
early-type stars where the disc density and tilt angle is highest. There is no
clear pattern to why those few models show stronger V/R ratios than others, so
it would be difficult to discern, without further constraints by other
observables, whether V/R ratios in actual observations of Be stars are due to
a tilted disc or a density enhancement in the disc, like those produced by
spiral enhancements in $\zeta$ Tau (Štefl et al., 2009; Escolano et al.,
2015), and 48 Lib (Silaj et al., 2016).

The second case is where the disc is already tilted, and undergoing precession
due to the influence of a misaligned binary companion. We are able to simulate
the disc precessing about the stellar pole by holding the $\theta$ observing
angle constant and changing $\phi$ only. Here we find that the percent
polarization, $V$ magnitude, and H$\rm\alpha$ EW oscillate at half the
precession period, although the oscillation will not be perfectly symmetric
unless the observer is directly aligned with the stellar equator. The position
angle, on the other hand, will oscillate in sync with disc precession. The V/R
ratio undergoes an antisymmetric half-period oscillation, with nodes at
$\phi\,=\,$$$ and $None$ due to the violet and red sides of the disc reversing
when the observer moves to the other side of the disc.

We then investigated two other scenarios. First, we computed two truncated
disc models out to a radius of 20 $R_{*}$, to see any possible obscuring
effects that may happen compared to a disc 50 $R_{*}$ in size. These
simulations revealed that the outer disc from 20 to 50 stellar radii only
marginally increased H$\rm\alpha$ emission, while not changing other examined
observables. The temperature structure was also unchanged out to 20 stellar
radii, as expected.

Second, and more importantly, we computed two models that were linearly warped
up to a maximum angle of $None$. This model revealed a cooler outer disc
temperature versus its flat tilted counterpart, and an inner temperature
structure that followed the warp of the disc. The observables of this model
are essentially a mix of all the tilted models together, with some observables
better matching the non-tilted or $None$ models, while others matched the
higher tilt models. This shows how important it is to recognize that Be star
discs emit some wavelengths from dense inner volumes while other wavelengths
come from larger radial positions in the disc. These warped models are merely
an initial test to see what effects a warped disc may introduce. A proper
warped disc study is beyond the scope of this paper, though certainly merits
its own study.

These simulations are a vital step towards simulations of more complex disc
configurations, such as ones containing warped discs, or those presented by
Suffak et al. (2022), which contain holes and tearing of the disc. The flat-
tilted models here will be a good benchmark for analysis of these discs that
are tilted, warped, and have asymmetric density distributions. With the
fundamentals presented here, we will be able to tackle more complicated Be
star systems such as Pleione, which is suspected to have a periodic tearing
disc (Marr et al., 2022; Martin & Lepp, 2022).

## Acknowledgements

We thank the anonymous referee for their thorough comments which improved the
paper. We greatly acknowledge the work of Marr (2022), whose preliminary work
on the temperature of tilted discs inspired and aided this work. C.E.J.
acknowledges support through the National Science and Engineering Research
Council of Canada. A. C. C. acknowledges support from CNPq (grant
311446/2019-1) and FAPESP (grants 2018/04055-8 and 2019/13354-1). THA
acknowledges support from FAPESP (grant 2021/01891-2 ). This work was made
possible through the use of the Shared Hierarchical Academic Research
Computing Network (SHARCNET).

## Data Availability

Although there is no observational data, the hdust models computed for this
work can be made available upon request.

## References

  * Baade et al. (2016) Baade D., et al., 2016, A&A, 588, A56
  * Brown et al. (2019) Brown R. O., Coe M. J., Ho W. C. G., Okazaki A. T., 2019, MNRAS, 488, 387
  * Carciofi & Bjorkman (2006) Carciofi A. C., Bjorkman J. E., 2006, ApJ, 639, 1081
  * Carciofi & Bjorkman (2008) Carciofi A. C., Bjorkman J. E., 2008, ApJ, 684, 1374
  * Collins (1987) Collins George W. I., 1987, in Slettebak A., Snow T. P., eds, IAU Colloq. 92: Physics of Be Stars. p. 3
  * Cox (2000) Cox A. N., 2000, Allen’s astrophysical quantities
  * Cyr et al. (2017) Cyr I. H., Jones C. E., Panoglou D., Carciofi A. C., Okazaki A. T., 2017, MNRAS, 471, 596
  * Escolano et al. (2015) Escolano C., Carciofi A. C., Okazaki A. T., Rivinius T., Baade D., Štefl S., 2015, A&A, 576, A112
  * Ghoreyshi et al. (2018) Ghoreyshi M. R., et al., 2018, MNRAS, 479, 2214
  * Ghoreyshi et al. (2021) Ghoreyshi M. R., Carciofi A. C., Jones C. E., Faes D. M., Baade D., Rivinius T., 2021, ApJ, 909, 149
  * Hanuschik et al. (1996) Hanuschik R. W., Hummel W., Sutorius E., Dietle O., Thimm G., 1996, A&AS, 116, 309
  * Jones et al. (2004) Jones C. E., Sigut T. A. A., Marlborough J. M., 2004, MNRAS, 352, 841
  * Jones et al. (2008) Jones C. E., Tycner C., Sigut T. A. A., Benson J. A., Hutter D. J., 2008, ApJ, 687, 598
  * Klement et al. (2015) Klement R., et al., 2015, A&A, 584, A85
  * Lee et al. (1991) Lee U., Osaki Y., Saio H., 1991, MNRAS, 250, 432
  * Marr (2022) Marr K., 2022, PhD thesis, University of Western Ontario, London, ON, CA, https://ir.lib.uwo.ca/etd/8376
  * Marr et al. (2021) Marr K. C., Jones C. E., Carciofi A. C., Rubio A. C., Mota B. C., Ghoreyshi M. R., Hatfield D. W., Rímulo L. R., 2021, ApJ, 912, 76
  * Marr et al. (2022) Marr K. C., Jones C. E., Tycner C., Carciofi A. C., Silva A. C. F., 2022, ApJ, 928, 145
  * Martin & Lepp (2022) Martin R. G., Lepp S., 2022, MNRAS, 516, L86
  * Martin et al. (2011) Martin R. G., Pringle J. E., Tout C. A., Lubow S. H., 2011, MNRAS, 416, 2827
  * Martin et al. (2014a) Martin R. G., Nixon C., Armitage P. J., Lubow S. H., Price D. J., 2014a, ApJ, 790, L34
  * Martin et al. (2014b) Martin R. G., Nixon C., Lubow S. H., Armitage P. J., Price D. J., Doğan S., King A., 2014b, ApJ, 792, L33
  * McGill et al. (2011) McGill M. A., Sigut T. A. A., Jones C. E., 2011, ApJ, 743, 111
  * McGill et al. (2013) McGill M. A., Sigut T. A. A., Jones C. E., 2013, ApJS, 204, 2
  * Millar & Marlborough (1998) Millar C. E., Marlborough J. M., 1998, ApJ, 494, 715
  * Millar & Marlborough (1999a) Millar C. E., Marlborough J. M., 1999a, ApJ, 516, 276
  * Millar & Marlborough (1999b) Millar C. E., Marlborough J. M., 1999b, ApJ, 516, 280
  * Millar & Marlborough (1999c) Millar C. E., Marlborough J. M., 1999c, ApJ, 526, 400
  * Millar et al. (2000) Millar C. E., Sigut T. A. A., Marlborough J. M., 2000, MNRAS, 312, 465
  * Rímulo et al. (2018) Rímulo L. R., et al., 2018, MNRAS, 476, 3555
  * Rivinius et al. (2013) Rivinius T., Carciofi A. C., Martayan C., 2013, A&ARv, 21, 69
  * Sigut & Jones (2007) Sigut T. A. A., Jones C. E., 2007, ApJ, 668, 481
  * Silaj et al. (2010) Silaj J., Jones C. E., Tycner C., Sigut T. A. A., Smith A. D., 2010, ApJS, 187, 228
  * Silaj et al. (2016) Silaj J., et al., 2016, ApJ, 826, 81
  * Suffak et al. (2020) Suffak M. W., Jones C. E., Tycner C., Henry G. W., Carciofi A. C., Mota B. C., Rubio A. C., 2020, ApJ, 890, 86
  * Suffak et al. (2022) Suffak M., Jones C. E., Carciofi A. C., 2022, MNRAS, 509, 931
  * Vieira et al. (2017) Vieira R. G., Carciofi A. C., Bjorkman J. E., Rivinius T., Baade D., Rímulo L. R., 2017, MNRAS, 464, 3071
  * Waters (1986) Waters L. B. F. M., 1986, A&A, 162, 121
  * Štefl et al. (2009) Štefl S., et al., 2009, A&A, 504, 929

##  Appendix A Tilted Temperature cross sections

Figure 17: Same as Figure 6, but for models 2-7.

Figure 18: Same as Figure 6, but for models 8-13.

Figure 19: Same as Figure 6, but for models 14-19.

Figure 20: Same as Figure 6, but for models 20-25.

Figure 21: Same as Figure 6, but for models 26-29, 31 and 32

##  Appendix B Midplane Temperatures of Tilted vs Non-Tilted Discs

Figure 22: Temperature vs. radius at the disc midplane for all models with a
B0 or B2 central star, in the direction $\phi\,=\,$$$. The four lines are for
four different tilt angles as indicated by the legend.

Figure 23: Same as Figure 22 but for our B5 and B8 models.