# Introducing two improved methods for approximating radiative cooling in
hydrodynamical simulations of accretion discs

Alison K. Young,1,2 Maggie Celeste,3 Richard A. Booth,4 Ken Rice,1,2 Adam
Koval,1,2 Ethan Carter,5 and Dimitris Stamatellos5
1SUPA, Institute for Astronomy, University of Edinburgh, Blackford Hill,
Edinburgh EH9 3HJ, UK
2Centre for Exoplanet Science, University of Edinburgh, Edinburgh, EH9 3HJ, UK
3Institute of Astronomy, University of Cambridge, Madingley Rd, Cambridge CB3
0HA
4 School of Physics and Astronomy, Sir William Henry Bragg Building, Woodhouse
Ln., University of Leeds, Leeds LS2 9JT, UK
5 Jeremiah Horrocks Institute for Mathematics, Physics, and Astronomy,
University of Central Lancashire, Preston PR1 2HE, UK
E-mail<EMAIL_ADDRESS>(AKY)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The evolution of many astrophysical systems depends strongly on the balance
between heating and cooling, in particular star formation in giant molecular
clouds and the evolution of young protostellar systems. Protostellar discs are
susceptible to the gravitational instability, which can play a key role in
their evolution and in planet formation. The strength of the instability
depends on the rate at which the system loses thermal energy. To study the
evolution of these systems, we require radiative cooling approximations
because full radiative transfer is generally too expensive to be coupled to
hydrodynamical models. Here we present two new approximate methods for
computing radiative cooling that make use of the polytropic cooling
approximation. This approach invokes the assumption that each parcel of gas is
located within a spherical pseudo-cloud which can then be used to approximate
the optical depth. The first method combines the methods introduced by
Stamatellos et al. and Lombardi et al. to overcome the limitations of each
method at low and high optical depths respectively. The second, the “Modified
Lombardi” method, is specifically tailored for self-gravitating discs. This
modifies the scale height estimate from the method of Lombardi et al. using
the analytical scale height for a self-gravitating disc. We show that the
Modified Lombardi method provides an excellent approximation for the column
density in a fragmenting disc, a regime in which the existing methods fail to
recover the clumps and spiral structures. We therefore recommend this improved
radiative cooling method for more realistic simulations of self-gravitating
discs.

###### keywords:

hydrodynamics – radiative transfer – methods:numerical

††pubyear: 2024††pagerange: Introducing two improved methods for approximating
radiative cooling in hydrodynamical simulations of accretion discs–Introducing
two improved methods for approximating radiative cooling in hydrodynamical
simulations of accretion discs

##  1 Introduction

Radiative transfer is an important aspect of many hydrodynamical simulations.
Accurate treatment of heating and cooling is particularly important for
systems undergoing gravitational instability or fragmentation, such as star-
forming clouds and young, self-gravitating protoplanetary discs, in which the
strength of the gravitational instability is a function of the thermal balance
between shock heating and radiative cooling (Gammie, 2001). This determines
whether the disc develops spiral arms, fragments or maintains a smooth,
axisymmetric structure.

Early work using isothermal simulations suggested that a self-gravitating
phase should be very short-lived, with the disc either fragmenting quickly or
spiral instabilities acting to stabilize the disc’s density profile (Laughlin
& Bodenheimer, 1994). However, simulations that implemented more detailed
thermodynamics demonstrated that this was mostly a consequence of the
isothermal assumption and that a self-gravitating phase may persist for many
orbits (Lodato & Rice, 2004). Now that we expect planet formation to begin
very early when discs are likely to be self-gravitating (Nixon et al., 2018),
we need robust models of self-gravitating discs to address key questions
including how planetesimals grow in young discs. However, directly
implementing radiative transfer in 3-D simulations is by no means trivial. A
full treatment would require wavelength-dependent and potentially anisotropic
absorption, scattering, and emission of photons; tracing of the photons’ paths
through the system; and variable opacities and emissivities for both dust and
gas components. Additionally, explicit radiative transfer timesteps are
prohibitively short relative to hydrodynamic timescales, particularly in
optically thick locations, greatly increasing computation time (Castor, 2004;
Brandenburg & Das, 2020). These difficulties have prompted the development of
a number of approximate radiative transfer methods, each depending on some
simplifying assumptions and approximations. Such methods are naturally limited
in accuracy and scope compared to a full radiative transfer scheme, but allow
for the simultaneous simulation of both hydrodynamic and radiative processes.

One such method which has been used to study the effects of cooling in
protoplanetary discs is the so-called $\beta$-cooling model, developed by
Gammie (2001). This approach parameterises the local cooling time $t_{c}$ in
terms of a $\beta$ parameter and the local Keplerian angular velocity
$\omega_{\rm K}$, such that $t_{c}=\beta/\omega_{\rm K}$. The simplicity of
this parameterisation and ease of implementation in hydrodynamic codes has
made it a popular choice for studying the relationship between cooling and
fragmentation in discs (e.g. Rice et al., 2003; Cossins et al., 2009; Meru &
Bate, 2011; Boss, 2017), but it does not paint a complete picture. There is no
compelling argument for the $\beta$ parameter remaining constant throughout
the extent of the disc (Mercer et al., 2018), and more realistic cooling
methods yield significantly different results (Johnson & Gammie, 2003;
Vorobyov et al., 2020).

There are two main classes of more physically-motivated cooling methods for
protoplanetary disc simulations: those that impose an external radiation field
and calculate radiation transport through the disc; and those that make use of
local parameters to estimate the thermal balance at a given location in the
disc. Examples of the first approach include codes that make use of the flux-
limited diffusion approximation (Levermore & Pomraning, 1981), such as Ayliffe
& Bate (2011) and Meru & Bate (2010) with fixed prescribed temperatures at the
edge of the model; or Monte-Carlo methods to transport photon “packets”
throughout the disc, such as the coupled mcfost + phantom model (Pinte et al.,
2006; Price et al., 2018; Nealon et al., 2020). For the coupled radiative
transfer and hydrodynamics approach, a radiative equilibrium calculation is
performed every few hydrodynamical timesteps (this is chosen to be as
infrequent as possible to reduce the overheads), after which the temperatures
are updated. Although the radiative transfer is more accurate than under the
flux-limited diffusion approximation and can, for example, reproduce the
effects of shadowing, the coupled method is not time-dependent due to the
radiative equilibrium assumption. The outcome of this is that temperatures are
likely to be under- or over-estimated in different regions. The second
approach – estimating the cooling rate based on local properties – is the
focus of this paper. Existing methods are sufficient for modelling symmetrical
systems, i.e. spherical cloud collapse, but are less accurate for modelling
discs since they overestimate the optical depth in the disc. This then changes
the cooling time, which in turn affects the evolution of self-gravitating
discs. The use of Monte Carlo radiative transfer would allow for irregular
structures but is very costly and currently imposes radiative equilibrium. As
discussed above, it is crucial to have a robust cooling model for studying
gravito-turbulent discs, therefore a more accurate approximation for the
optical depth within discs and therefore the cooling rate is needed, without
adding a large computational expense.

In section 2 we first describe the methods that are currently used to
approximate radiative cooling by assuming that each parcel of gas, or SPH
particle, is embedded within a spherical gas cloud (“polytropic cooling”) and
using this to estimate the optical depth. The limitations of these methods are
detailed in 2.5 and we further motivate the development of our new techniques.
In section 3.1 we introduce the first of our improved methods for
approximating radiative cooling in discs, which combines the two. We present
our second method in section 3.2 which provides an alternative means of
regularising the Lombardi et al. (2015) method in the mid-plane by making use
of disc-specific geometry.

##  2 Existing polytropic cooling methods

Much of what we present here will be based on simulations using Smoothed
Particle Hydrodyamics (SPH) (e.g., Monaghan, 1992), but it is also applicable
to other hydrodynamic formalisms. In SPH, each particle represents a parcel of
gas and to estimate its radiative cooling rate we need to know its optical
depth. Physically, this involves integrating the opacity and density along the
path from the particle to the edge of the particle distribution (e.g., in the
case of a disc, to the disc surface). The cooling approximations described
below instead estimate a mean column density from the particle to the surface
and combine this with mean opacity estimates. Two different approaches for
estimating this mean column density are described, firstly using the
gravitational potential (Stamatellos et al., 2007) and secondly using the
pressure gradient (Lombardi et al., 2015).

###  2.1 The Stamatellos Method: The gravitational potential approach

Stamatellos et al. (2007) introduced a method for approximating the optical
depth of an SPH particle from its gravitational potential. The gravitational
potential energy is already calculated in SPH simulations that include the gas
self-gravity, hence minimising additional computational expense. This approach
treats each parcel of gas as if it were embedded within a spherically
symmetric ‘pseudo-cloud’ with polytropic density and temperature profiles. The
scale-length of the polytrope is determined by the gravitational potential at
the location of the SPH particle. Because we don’t know where this gas parcel
lies within the pseudo-cloud, we must calculate a mass-weighted average of the
pseudo-column density over all possible positions within it. The mean column
density for particle $i$ is then given by

$\bar{\Sigma}_{i}=\zeta_{n}\left[\frac{-\psi_{i}\rho_{i}}{4\pi
G}\right]^{1/2},$ (1)

where $\zeta_{n}$ is a constant that depends on the polytropic index (which is
related to $n$; see Stamatellos et al. 2007 for the full derivation), and
$\psi_{i}$ and $\rho_{i}$ are the gravitational potential and density at the
position of the particle respectively. Following Stamatellos et al. (2007) and
Forgan et al. (2009) we use $\zeta_{2}=0.368$. We now require a pseudo-mean
mass opacity $\bar{\kappa_{i}}(\rho,T)$, to estimate the mass weighted optical
depth. The pseudo-mean opacity for particle $i$ is calculated by averaging
over locations in the pseudo-cloud in the same manner as for the mean column
density. This allows us to estimate the optical depth using

$\bar{\tau_{i}}=\bar{\Sigma}_{i}\bar{\kappa_{i}}.$ (2)

###  2.2 The Lombardi Method: the pressure gradient approach

A slightly different approach was introduced by Lombardi et al. (2015) in
which the pressure scale-height is used to estimate the optical depth rather
than the gravitational potential. This method still assumes that each particle
is embedded within a polytropic pseudo-cloud but by using the pressure
gradient at the location of the particle, it performs better for non-spherical
gas distributions.

The pressure scale-height, $H_{{\rm p},i}$, of particle $i$ is given by:

$H_{{\rm p},i}\equiv P_{i}/\left|\nabla P_{i}\right|,$ (3)

where $P_{i}$ is the pressure at the location of particle $i$, and $\nabla
P_{i}$ is the pressure gradient.

The hydrodynamic acceleration excluding the contributions from artificial (and
real) viscosity and gravity is

$\boldsymbol{a}_{{\mathrm{h}},i}=-\frac{\nabla P_{i}}{\rho_{i}}.$ (4)

In smoothed particle hydrodynamics, the hydrodynamic acceleration
$\boldsymbol{a}_{{\mathrm{h}},i}$ is already computed during the force
calculation. For particle $i$, this is calculated via the following sum over
$j$ neighbour particles:

$\boldsymbol{a}_{{\mathrm{h}},i}\equiv\frac{d\boldsymbol{v}_{i}}{dt}=-\sum_{j}\left[\frac{m_{j}P_{i}}{\rho_{i}^{2}\Omega_{i}}\nabla
W_{i}+\frac{m_{j}P_{j}}{\rho_{j}^{2}\Omega_{j}}\nabla W_{j}\right],$ (5)

in which $m_{j}$ is the particle mass, $\rho_{i}$ ($\rho_{j})$ is the gas
density evaluated at particle $i$ ($j$), and the gradient of the smoothing
kernel is $\nabla W$. $\Omega$ is related to the gradient of the smoothing
length $h_{i}$

$\Omega_{i}\equiv 1-\frac{\partial
h_{i}}{\partial\rho_{i}}\sum_{j}m_{j}\frac{\partial W_{ij}(h_{i})}{\partial
h_{i}}.$ (6)

We can then obtain the pressure scale-height via

$H_{{\rm
p},i}=\frac{P_{i}}{\rho_{i}\left|\boldsymbol{a}_{{\mathrm{h}},i}\right|}.$ (7)

Calculating the mass-weighted pseudo-mean column density in the same manner as
Stamatellos et al. (2007), Lombardi et al. (2015) obtain

$\bar{\Sigma}_{i}=\zeta^{\prime}\rho_{i}H_{{\rm
P},i}=\frac{\zeta^{\prime}P_{i}}{\left|\boldsymbol{a}_{{\rm h},i}\right|},$
(8)

with $\zeta^{\prime}=1.014$. The optical depth can then be estimated using
Equation 2.

###  2.3 Equation of state

We employ an identical approach to Stamatellos et al. (2007) in which values
of the pseudo-mean opacity $\bar{\kappa_{i}}$, Planck mean opacity
$\kappa_{i}$, mean molecular weight $\mu_{i}$, and specific internal energy
$u_{i}$ are pre-computed for the required density and temperature ranges and
stored in a look-up table to reduce computational expense. These were
calculated following Bell & Lin (1994). The Rosseland-mean and Planck-mean
opacities are assumed to be the same and to depend only on density and
temperature and are given by:

$\kappa_{\rm R}(\rho,T)=\kappa_{\rm P}(\rho,T)=\kappa_{0}\rho^{a}T^{b},$ (9)

where the constants $\kappa_{0}$, $a$ and $b$ are chosen for each main
physical process contributing to the total opacity in that temperature and
density range (see Bell & Lin 1994; Stamatellos et al. 2007 for further
details). The tabulated values are then interpolated during the live
hydrodynamical computation and the ideal gas approximation is used to obtain
the gas pressure.

###  2.4 Updating the energy

Both methods give us a mass-weighted column density and tabulated values of
the pseudo-mean opacity. These quantities are now used to calculate the
radiative heating rate of gas parcel represented by particle $i$:

$\frac{du_{i}}{dt}\Big{|}_{\rm
rad}=\frac{4\sigma\left(T_{0}^{4}(\boldsymbol{r}_{i})-T_{i}^{4}\right)}{\bar{\Sigma}^{2}\bar{\kappa_{i}}(\rho_{i},T_{i})+\kappa_{i}^{-1}(\rho_{i},T_{i})},$
(10)

where $T_{0}$ is a minimum temperature determined by background radiation and
$\kappa_{i}^{-1}(\rho_{i},T_{i})$ is the Planck mean opacity, which is
tabulated along with the mass-weighted opacity
$\bar{\kappa_{i}}(\rho_{i},T_{i})$. External heating from nearby stars can be
included by calculating $T_{0}$ for each gas particle to prevent the gas
cooling radiatively below this minimum temperature.111Lombardi et al. (2015)
introduce a slightly different method for calculating the radiative heating
rate but we opt to use the original approach of Stamatellos et al. (2007) for
both methods.

The equilibrium temperature of particle $i$ is found by assuming that the
radiative cooling and the compressive and viscous heating are in balance:

$\frac{du_{i}}{dt}\Big{|}_{\rm hydro}+\frac{du_{i}}{dt}\Big{|}_{\rm rad}=0.$
(11)

Substituting the equilibrium temperature $T_{{\rm eq},i}$ for $T_{i}$ in
Equation 10 gives:

$T_{\rm{eq},i}^{4}=\frac{1}{4\sigma}\left(\bar{\Sigma}^{2}\bar{\kappa_{i}}(\rho_{i},T_{i})+\kappa_{i}^{-1}(\rho_{i},T_{i})\right)\frac{du_{i}}{dt}\Big{|}_{\rm
hydro}+T_{0}^{4}(\boldsymbol{r}_{i}).$ (12)

In the case that $T_{\rm{eq},i}^{4}\leq T_{0}^{4}(\boldsymbol{r}_{i})$, we set
$T_{\rm{eq},i}^{4}=T_{0}^{4}(\boldsymbol{r}_{i})$. The equilibrium energy is
then $u_{{\rm eq},i}=u(T_{{\rm eq,}i},\rho_{i})$. The energy is updated by
considering the approach towards the equilibrium temperature with thermal
timescale

$t_{{\rm
therm},i}=\left(u_{\rm{eq},i}-u_{i}\right)\left[\frac{du_{i}}{dt}\Big{|}_{\rm
hydro}+\frac{du_{i}}{dt}\Big{|}_{\rm rad}\right]^{-1}.$ (13)

For implementation within a typical SPH code, the cooling rate is then given
by

$\frac{du_{i}}{dt}\Big{|}_{\rm cool}=\frac{1}{\delta
t}\left[u_{i}\exp\left(\frac{-\delta
t}{t_{\rm{therm},i}}\right)+u_{\rm{eq},i}\left(1-\exp\left(\frac{-\delta
t}{t_{\rm{therm},i}}\right)\right)\right]\\\ +\frac{du_{i}}{dt}\Big{|}_{\rm
hydro}.$ (14)

The energy equation is then updated via

$u_{i}(t+\delta t)=u_{i}(t)+\delta t\frac{du_{i}}{dt}\Big{|}_{\rm cool}.$ (15)

###  2.5 Limitations of these methods

In both cases, given the estimate is based only on local parameters already
used in hydrodynamic calculations in SPH, there is no need for communication
of information with other fluid elements, and there is therefore a minor
computational cost when compared with methods that utilise a radiation field
or photon packets. However, any local approximation of optical depth will
necessarily require assumptions about the system’s global geometry. In the
case where one wishes to study systems with multiple geometries co-occurring –
for example, discs with spherical clumps and spiral waves, or with low-density
envelopes or streamer tails – accurately estimating the cooling rate in all
regions of the disc with just one method can be problematic.

The Stamatellos et al. (2007) method significantly over-estimates the optical
depth in lower density regions of the disc but performs better towards the
disc mid-plane and in spherical clumps. The Lombardi et al. (2015) method is
more accurate overall in disc simulations (Mercer et al., 2018), but breaks
down in regions where the pressure gradient approaches zero, such as in the
disc mid-plane or the center of clumps. By combining the two methods, one can
take advantage of the strengths of each approach to achieve more accurate
cooling throughout the whole of the disc. It is also possible to improve on
the Lombardi et al. (2015) method in the regions where it fails by using
another means of estimating the optical depth that makes use of a disc-
specific density distribution assumption, as there is a simple relationship
between mid-plane density, scale-height and column density in discs.

##  3 Method

In this section we detail two new methods for estimating the optical depth of
a particle. First, in 3.1 we combine the estimates from the methods of
Stamatellos et al. (2007) and Lombardi et al. (2015), and second in 3.2 we
modify the method of Lombardi et al. (2015) specifically tailored for
modelling self-gravitating discs.

###  3.1 The Combined Method: combining estimates of optical depth

Figure 1: The column density from $R$ to $R_{\rm out}$ for an $n=3/2$
polytrope (blue line, "Exact") compared with the estimate from Lombardi et al.
(2015). Also shown are the estimate generated by the combined method (see
subsection 3.1) and a modification of the Lombardi method for discs (see
subsection 3.2). Figure 2: Top: The column density from $z$ to infinity for
the analytic self-gravitating disc model. Also shown are different estimates
based on the Lombardi method and Combined method, as well as a Modified
Lombardi method applicable to discs (see subsection 3.2). Bottom: The ratio of
the estimated column density to the exact column density.

While Lombardi et al. (2015) and Mercer et al. (2018) showed that the pressure
scale-height and density can be used to provide a good estimate of the column
density (Equation 8), this estimate breaks down in regions where the pressure
gradient is close to zero, such as at the centre of a fragment or close to a
disc mid-plane. Since a large part of the mass can be found under such
conditions, the average cooling rate can be affected and we are motivated to
find a way to improve the estimate under these conditions.

Here we explore replacing the Lombardi method with a surface density estimated
by the Stamatellos method in regions where $\nabla P\rightarrow 0$ as the
Stamatellos method remains well defined in such regions.

First, we define a scale-height based on the Stamatellos method via

$H_{\rm S}=\frac{\zeta_{n}}{\zeta^{\prime}}\left[\frac{-\psi}{4\pi
G\rho}\right]^{1/2},$ (16)

which we combine with the Lombardi estimate ${H_{\rm P}}$ in inverse
quadrature, i.e.

$H_{\rm C}=\left(H_{\rm P}^{-2}+H_{\rm S}^{-2}\right)^{-1/2}.$ (17)

This selects the smaller of the two estimates, chosen because the Stamatellos
method is known to overestimate the surface density in disc geometries and
therefore generally we will only want to use this estimate when the Lombardi
method is also failing, producing overestimates itself.

As a first test of this Combined method we consider a polytrope as an
approximation to the structure of a self-gravitating fragment, and show
$\Sigma(z)$ for an $n=3/2$ polytrope and also the different column density
estimates in Figure 1. As previously stated, the Lombardi method works well
away from the centre of the polytrope but overestimates the column density
once the density and pressure gradients tend to zero close to the centre of
the polytrope. The Combined method offers an improvement over the Lombardi
method alone near the core, but it does lead to a moderate underestimate of
the column. This follows because the Stamatellos method uses an estimate of
the column density that is averaged over the entire polytrope and therefore
must underestimate the column at the centre, where the column is highest.

Next, we consider how the Combined method works in a disc geometry. To this
end, we consider a thin, self-gravitating disc. High-resolution simulations
show that self-gravitating discs have a temperature structure that is
typically close to vertically isothermal (Shi & Chiang, 2014; Booth & Clarke,
2019), so we choose a vertically-isothermal structure. Then vertical
hydrostatic equilibrium implies

$\frac{1}{P}\frac{{\rm d}P}{{\rm d}z}=\frac{1}{\rho}\frac{{\rm d}\rho}{{\rm
d}z}=-\frac{\omega_{\rm K}^{2}z+4\pi G\rho_{0}\tilde{\Sigma}}{c_{\rm
s}^{2}}\\\ =-\frac{1}{H_{*}^{2}}\left(z+\frac{\tilde{\Sigma}}{Q_{\rm
3D}}\right)=-H_{P}^{-1},$ (18)

where we have assumed that the disc’s self-gravity can be calculated using a
slab approximation, as in Wilkins & Clarke (2012) for example. Here
$H_{*}=c_{\rm s}/\omega_{\rm K}$, where $c_{\rm s}$ is the isothermal sound
speed, $\omega_{\rm K}$ is the Keplerian angular frequency, and
$\tilde{\Sigma}=\int_{0}^{z}\rho(z)/\rho(0){\rm d}z$. For the last equality we
have introduced $Q_{\rm 3D}=\omega_{\rm K}^{2}/4\pi G\rho(0)$, which is the 3D
definition of the Toomre Q parameter (e.g. Mamatsashvili & Rice, 2010; Lin &
Kratter, 2016).

To compute the potential we need to extend this 1D model to a global disc
because the potential of an infinitely extended self-gravitating slab is
undefined. We choose a disc that has a constant $Q_{\rm 3D}$ with radius so
the vertical structure is given by the same solution to Equation 18
everywhere. The definition of $Q_{\rm 3D}$ then enforces the mid-plane density
scaling as $\rho_{0}\propto R^{-3}$. Specifying the scale-height, $H_{*}$,
fixes the surface density and the disc model is complete. Here we choose
$H_{*}/R=0.1$. We then compute the potential by direct integration of Green’s
function for Poisson’s equation in 3D. We assume that the disc extends from 5
au to 100 au for consistency with the discs used in the SPH simulations. We
note that while the final potential is not perfectly consistent with the disc
structure, the difference is small when the disc is thin (e.g. Bertin &
Lodato, 1999; Wilkins & Clarke, 2012)

In Figure 2 we compare the surface density for the analytic disc model to the
Lombardi estimate using the exact $H_{\rm P}$ from Equation 18. The estimate
works well at large $z$ but, as expected, it overestimates the surface density
at $z=0$. To regularize the Lombardi method we evaluate the potential assuming
a point 50 au from the star. Here, we see that the Combined method improves on
the Lombardi estimate, as expected. In this case, the Combined method
moderately overestimates the column at the mid-plane, which is expected since
the Stamatellos method has been shown previously to overestimate the column in
disc geometries (Wilkins & Clarke, 2012; Lombardi et al., 2015; Mercer et al.,
2018).

These two simple tests demonstrate that combining the Lombardi and Stamatellos
methods should offer a significant improvement over using either method alone.
In section 4 we explore this for more realistic scenarios that arise in SPH
simulations.

###  3.2 The Modified Lombardi method: An improved Lombardi method for disc
simulations

The Stamatellos method is not the only viable option for regularizing the
Lombardi method. Since the Lombardi method only fails in regions where $\nabla
P\approx 0$, such as the mid-plane of discs or the centre of polytropes, we
only need to replace the estimate at such locations. Here we explore an
alternative that makes use of the fact that there is a simple relationship
between the mid-plane density, scale-height and column density in a disc
geometry.

For a non-self-gravitating disc, the (half) column density and mid-plane
density are related via $\Sigma/2=\rho_{0}\sqrt{\pi/2}H_{*}$, and so
$H_{0}=\sqrt{\pi/2}H_{*}$ could be used in the place of $H_{\rm S}$. Away from
the mid-plane $H_{\rm P}<\sqrt{\pi/2}H_{*}$, so the Lombardi method would take
over. However, using $\sqrt{\pi/2}H_{*}$ overestimates the column density for
self-gravitating discs or in fragments because self-gravity further compresses
the gas, reducing the scale-height. Fortunately, it is possible to estimate
$H_{0}$ for a disc with any $Q_{\rm 3D}$. For $Q_{\rm 3D}\gg 1$,
$H_{0}=\sqrt{\pi/2}H_{*}$ suffices. In the regime $Q_{\rm 3D}\ll 1$ the second
term in Equation 18 dominates and the equation can be made dimensionless by
substituting $z$ with $z^{\prime}\lambda$ where $\lambda=H_{*}\sqrt{Q_{\rm
3D}}$. We therefore expect $H_{0}/H_{*}$ to scale with $\sqrt{Q_{\rm 3D}}$
when $Q_{\rm 3D}$ is small, and indeed find the fitting formula,

$\frac{H_{0}}{H_{*}}=\frac{t}{\sqrt{1+1/(tQ_{\rm 3D})}},$ (19)

describes $\Sigma(0)/\rho(0)$ well, where $t=\sqrt{\pi/2}$ (Figure 3) 222In
the SPH model, we implement $\omega_{{\rm
K},i}=\sqrt{GM_{*}/|\boldsymbol{r}_{K,i}|^{3}}$ where
$\boldsymbol{r}_{K,i}=\boldsymbol{r}_{i}-\boldsymbol{r}_{*}$ is the separation
between the gas particle $i$ and the sink particle (star). If there is no sink
particle, $\omega_{{\rm K},i}=0$. This approach would need to be adapted for
simulations with multiple stars, for example by adding $\omega_{{\rm K},i}$
for each sink particle in quadrature..

An advantage of Equation 19 is that it can be evaluated using quantities local
to the disc mid-plane only (i.e., $\rho$, $c_{\rm s}$, and $\omega_{\rm K}$).
Further, while Equation 19 is only strictly meaningful at the mid-plane it can
be evaluated anywhere as it only depends on local quantities. Since $H_{0}$
increases away from the mid-plane (because $\rho$ decreases) and $H_{\rm P}$
decreases, replacing $H_{\rm S}$ with Equation 19 in Equation 17 results in a
method that produces similar estimates to the Lombardi method away from the
mid-plane, which is generally quite accurate. Indeed, Figure 2 shows that this
‘Modified Lombardi’ method outperforms the normal Lombardi and Combined
Lombardi + Stamatellos methods in discs. For the method to be robust, however,
we need it to work well when the geometry deviates significantly from a disc.

Next, we apply this Modified Lombardi method to a polytrope as a model for a
fragment in a self-gravitating disc. Under such conditions, we expect $Q_{\rm
3D}\ll 1$, such that $H_{0}\approx 2c_{\rm s}/\sqrt{4\pi G\rho_{0}}$. We apply
this estimate of $H_{0}$ to the polytrope in Figure 1, which shows that
Equation 19 provides an estimate of $\Sigma(z)$ that is accurate to within a
factor of 2, in a spherical geometry. Since the size of a polytrope scales
with $c_{\rm s}\sqrt{(n+1)/4\pi G\rho(0)}$ (e.g. Clarke & Carswell, 2014), the
accuracy of this estimate is independent of $Q_{\rm 3D}$ as long as $Q_{\rm
3D}\ll 1$. Since the modified Lombardi method works well for disc geometries
in general and for polytropic clumps when $Q_{\rm 3D}\ll 1$, we expect it to
be widely applicable to the structures present in self-gravitating discs and
it also works well for spherical geometries. Finally, it should be noted that
in deriving Equation 19 we have assumed hydrostatic equilibrium. As a result,
the Modified Lombardi method may not offer a significant improvement away from
hydrostatic equilibrium.

Figure 3: A comparison between the half thickness of a self-gravitating disc
for different $Q_{\rm 3D}$ and the fit formula given in Equation 19.

##  4 Tests

We now present tests using the SPH code phantom (Price et al., 2018). For
comparisons with previous work, we first simulate the collapse of a uniform
density sphere and then consider both low- and high-mass discs. For the
collapse of the uniform density sphere, we only consider the Combined method,
but for the low- and high-mass protoplanetary discs we consider both the
Combined method and the Modified Lombardi method. The key factor that
determines how accurately the cooling rate is estimated and the difference
between the methods is how well the column density at each position is
estimated. The column density estimates are therefore the focus of this
section.

###  4.1 Collapse of a uniform density sphere

The gravitational collapse of a uniform density sphere represents the collapse
of a protostellar cloud core and forms a useful benchmark (see Masunaga &
Inutsuka 2000 for the first radiation hydrodynamic simulation of such a
collapse). For this test, we use $2\times 10^{6}$ SPH particles to simulate a
cloud core of mass $1.5$ M⊙, with an initial radius of $10^{4}$ au, and an
initially uniform temperature of $5$ K.

Fig. 4 shows the evolution of the maximum density and temperature of the cloud
core using the Stamatellos method, the Lombardi method, and the Combined
method. For the Stamatellos and Combined methods, this is simply the density
and temperature of the SPH particle with the highest density. For the Lombardi
method there is considerable scatter so we average the 200 SPH particles with
the highest density. The evolution of models with the Stamatellos, Lombardi
and Combined methods are very similar. Moreover, Figure 4 shows that these
methods produce results that are close to what is expected (black dash-dot
line from Masunaga & Inutsuka 2000). As mentioned in Stamatellos et al.
(2007), the differences with respect to Masunaga & Inutsuka (2000) are likely
due to the different opacities implemented.

Figure 4: The evolution of the maximum density and temperature during the
collapse of a uniform density sphere with three radiative cooling methods:
Stamatellos method (yellow line), Lombardi method (red dotted line), and the
Combined method (blue dashed line). The results from the simulation of
Masunaga & Inutsuka (2000) are plotted (black dashed-dot line) for comparison.
The initial mass and radius are $1.5$ M⊙ and $10^{4}$ au respectively, and the
temperature is initially 5 K throughout. The evolution under the Combined
method lies almost exactly on top of the results obtained using the
Stamatellos method.

The small differences in the evolution between the models are due to the
methods used to estimate the column density. Therefore, we now examine the
column density estimates in each case. For this analysis we used a single
simulation (using the Combined cooling method) and determined the column
density estimate for all three cooling methods at different times in this
simulation.

Figure 5: Estimates of the column density for a collapsing sphere with an
initially uniform density at $t=80000$ years (left-hand panel) and at
$t=145000$ years (right-hand panel). Estimates are calculated using the
Stamatellos method (yellow open circles), the Lombardi method (red open
squares), the Combined method (blue stars), and the Modified Lombardi method
(green triangles). The dashed line is the “true” column density profile,
determined by integrating the density in the SPH simulation from the centre of
the cloud to the surface.

Figure 5 shows the resulting column density estimates as a function of radius
for two snapshots during the collapse. We compare these estimates with the
“true” column density (black dashed lines), calculated by integrating the
density outwards from the specified cloud radius to the surface. Near the
cloud core centre, the Lombardi method initially over-estimates the column
density and Fig. 5 displays a lot of scatter. This is because the pressure
gradient will initially be close to zero near the core centre and so the
Lombardi method will not produce reliable estimates of the column density. It
is because of this scatter that it is necessary to average the 200 particles
with the highest density when producing the results for the Lombardi shown in
Fig. 4.

On the other hand, near the cloud centre, the Stamatellos method, the Combined
method (which tends to the Stamatellos method in these regions) and the
Modified Lombardi method initially under-estimate the column density, but
avoid the large scatter seen in the Lombardi method. However, at larger radii,
the Stamatellos method clearly over-estimates the column density, while the
Lombardi and Modified Lombardi methods produce values much closer to the real
column density. In these regions, the Combined method tends to the Lombardi
method and produces a reliable estimate of the column density.

The right-hand panel of Figure 5 shows that as the cloud collapse progresses,
the Combined method tends to the Lombardi method throughout most of the cloud
core and both return reasonable estimates of the column density.The Modified
Lombardi method also provides a reasonable estimate of the column density,
despite being designed for a disc geometry.

###  4.2 Evolution of a protoplanetary disc

####  4.2.1 Low-mass disc

We first consider a low-mass disc, similar to that presented in Mercer et al.
(2018). The disc has an initial surface density profile of $\Sigma\propto
r^{-1}$, extending from $r=5$ au to $r=100$ au, and an initial sound speed
profile of $c_{\rm s}\propto r^{-0.25}$. The total disc mass is $M_{\rm
disc}=0.01$ M⊙, and the aspect ratio at $r=5$ au is initially $H/r=c_{\rm
s}/\Omega=0.05$. The disc gas is represented by $2\times 10^{6}$ SPH particles
and the central star has a mass of $M_{*}=0.8$ M⊙ and is represented by a sink
particle. The artificial viscosity parameters are $\alpha_{\rm SPH}$ = 1.0 and
$\beta_{\rm SPH}=2.0$.

We run the simulation for $t=1200$ years to allow the disc to settle and
specify a minimum disc temperature of $T_{0}=25$ K. The column density
snapshot of the disc at $t=1200$ years (Fig. 6) shows the disc has a smooth
structure. We then consider a thin azimuthal ring within the disc at $r=75$
au. Fig. 7 shows the vertical density structure at $r=75$ au with a by-eye fit
to the particle plot. Fitting with the expected profile of
$\rho=\rho_{o}\exp(-z^{2}/2H^{2})$ demonstrates that the scale-height at
$r=75$ au is about $H=6.7$ au.

Figure 6: Snapshot of the disc column density at $t=1200$ years for the
$M_{\rm disc}=0.01$ M⊙ (low-mass) disc discussed in Section 4.2.1. Figure 7:
By-eye fit to the vertical profile, at $r=75$ au, of the SPH simulation of a
0.01 M⊙ (low-mass) disc around a 0.8 M⊙ central star.

We now compare this value obtained directly from the simulation to estimates
obtained using each of the Stamatellos, Lombardi, Combined, and Modified
Lombardi methods. Estimates of the scale-height, $H$, for a randomly selected
sample of SPH particles are presented in Fig. 8, along with the value obtained
directly from the simulation ($H=6.7$ AU, dotted line, c.f., Fig. 7). Near the
midplane, the Stamatellos method (open yellow circles) produces a reasonable
estimate for the scale-height, but at higher altitudes is clearly far too
large. The Lombardi method is, again, unreliable near the midplane, since this
is a region where the pressure gradient is tending towards zero, but appears
to produce more reasonable estimates at higher altitudes. Combining the two
methods (Fig. 8, blue stars) then produces estimates that tend to the
Stamatellos method near the midplane, and to the Lombardi method at higher
altitudes. The Modified Lombardi method is an almost exact match to the
expected scale-height in the midplane, and then follows the Lombardi method at
higher altitudes.

Figure 8: Estimates for the disc scale-height at $r=75$ AU of the low-mass
disc calculated with the Stamatellos method (yellow circles), the Lombardi
method (red open squares), the Combined method (blue stars) and the Modified
Lombardi method (green triangles).The dotted line is $H=6.7$ AU, the “true”
scale-height obtained from the fit shown in Figure 7.

Fig. 9 shows the column density estimates found using equations 1 and 8 at
$r=75$ AU as a function of the height above the midplane, $z$, for the four
methods together with the full surface density from the SPH simulation (dotted
line) and column densities determined directly from the SPH simulation (dashed
line). Again, for the Lombardi, Stamatellos, Combined, and Modified Lombardi
methods, we plot a random sample of the estimates. The full SPH surface
density (dotted line) is determined by integrating the SPH density at $r=75$
AU over the full disc column, while the column density against $z$ (dashed
line) is determined by integrating the SPH density from $z$ to the nearest
disc surface.

The column density estimates from the Stamatellos, Lombardi, and Combined
methods are determined by multiplying the scale-height, $H$, in Figure 8 by
the local gas density. It should be noted that the column density estimates
from the four methods are one-sided column densities (i.e., from the $z$ to
the nearest disc surface) and, hence, at the midplane are about a factor of
two smaller than the full column/surface density (dotted line in Fig. 9).

Comparing the column densities from the SPH simulation itself with the various
estimates, shows that the Lombardi method is again unreliable near the
midplane ($z\sim 0$) but produces reliable estimates at higher $z$, while the
Stamatellos method produces reasonable estimates near the midplane, but
produces column density estimates that are much too high at higher $z$. Both
the Combined method and the Modified Lombardi method produce reasonable
estimates of the column density throughout the disc column.

Figure 9: Estimates for the column density against $z$ at $r=75$ AU of the
$0.01$ M⊙ (low-mass) disc using the Stamatellos method (yellow open circles),
the Lombardi method (red open squares), the Combined method (blue stars) and
the Modified Lombardi method (open green triangles). The dashed line shows the
column density determined directly from the simulation by integrating the
density from $z$ to the nearest disc surface. The total column density at
$r=75$ AU (equivalent to the surface density) is indicated by the dotted line
and will be twice the column density from the midplane to the nearest surface.
We are looking to obtain estimates as close to the dashed line as possible.

Finally, Fig. 10 shows the midplane column density estimates $\bar{\Sigma}(r)$
and half the disc surface density determined directly from the SPH simulation,
which we take as a reasonable estimate of the column density from the disc
midplane to the surface. For the Lombardi, Stamatellos, Combined and Modified
Lombardi method estimates, we plot a random sample of column density estimates
that are within $0.05$ AU of the midplane ($z=0$). Again, the Lombardi method
shows a lot of scatter, while the combined method removes most of this scatter
and produces reasonable estimates for the column density at almost all radii.
However, it is clear that the Stamatellos, Lombardi, and Combined methods all
slightly over-estimate the midplane column densities, with the discrepancy
becoming large at small disc radii. The Modified Lombardi method matches the
expected column density at all radii, illustrating how this is a particularly
reliable method for disc-like geometries.

Figure 10: Column density estimates against disc radius for the $0.01$ M⊙
(low-mass) disc. For the Stamatellos (yellow open circles), Lombardi (red open
squares), Combined (blue stars) and Modified Lombardi (open green triangles)
methods, we plot a random sample of estimates that are within $0.05$ AU of the
midplane ($z=0$).

####  4.2.2 High-mass disc

We now consider two different high-mass discs, both evolved using the Combined
method. Both discs initially extend from $r=5$ AU to $r=100$ AU, have surface
density profiles of $\Sigma\propto r^{-1}$, sound speed profiles of $c_{\rm
s}\propto r^{-1/4}$, initial scale-heights at $r=5$ AU of $H=0.275$ AU (i.e.,
$H/r=0.055$), and have central star masses of $M_{*}=0.8$ M⊙. The disc masses
are $M_{\rm disc}=0.2$ M⊙ and $M_{\rm disc}=0.25$ M⊙, and each disc is
represented by $2\times 10^{6}$ SPH particles. The artificial viscosity
parameters are $\alpha_{\rm SPH}$ = 1.0 and $\beta_{\rm SPH}=2.0$ and a floor
temperature of $T_{0}=10$ K was imposed. Fig. 11 shows that the $M_{\rm
disc}=0.2$ M⊙ settles into a quasi-steady state with spirals, while the
$M_{\rm disc}=0.25$ M⊙ disc becomes violently unstable and undergoes
fragmentation.

Figure 11: Column density renderings of the SPH simulations of high-mass discs
with masses of $M_{\rm disc}=0.2$ M⊙ (left-hand panel) and $M_{\rm disc}=0.25$
M⊙ (right-hand panel). Both were evolved using the Combined method.

We now consider how well the different methods for estimating the column
density represent structures within the disc. Fig. 12 shows the midplane
column density along a thin strip along the positive $x$-axis of the $M_{\rm
disc}=0.2$ M⊙ disc shown in the left-hand panel of Fig. 11. The region
considered is $x>0$ AU, $|y|<1$ AU. The black crosses show the “true” column
density determined by integrating the gas densities in the SPH simulation,
while the symbols are the estimates from the Lombardi, Stamatellos, Combined
and Modified Lombardi methods for $|z|<0.1$ AU.

Since the midplane is a region where the pressure gradient will tend to zero,
the Lombardi method again shows a lot of scatter. The Combined method (blue
stars) removes most of this scatter and seems to do reasonably well at
capturing the spiral structure in the disc, for example it captures the trough
at $r\sim 45$ AU and the peak at $r\sim 55$ AU. However, as in Fig. 10, the
Stamatellos, Lombardi and Combined methods tend to over-estimate the midplane
column density, generally by a factor of about 2. The Modified Lombardi method
not only captures the midplane density structures, but also mostly matches the
expected column density.

Figure 12: Column density against radius along a thin strip on the $x$-axis of
the $M_{\rm disc}=0.2$ M⊙ disc shown in the left-hand panel of Figure 11,
close to the midplane. The black crosses show the column density determined by
integrating the gas density in the SPH simulation, while the other symbols are
column densities from the Stamatellos method (yellow open circles), the
Lombardi method (red open squares), Combined method (blue stars), and Modified
Lombardi method (open green triangles).

Turning now to fragmenting discs, we examine how well the methods can recover
the structure of clumps within a disc. The clump considered here formed in the
$M_{\rm disc}=0.25$ M⊙ disc (Fig. 11, right panel) and is located at $x\sim
45$ AU, $y\sim 45$ AU in the snapshot. The column density as a function of the
radius from the centre of this clump is presented in Fig. 13. Since we don’t
know the exact boundary of the clump, we calculate the “true ” column density
(black dashed line) by integrating the gas density out to a clump radius of
$10$ AU, which could slightly over-estimate the column density as it may
include some of the surrounding gas disc. The Lombardi method over-estimates
the column density near the clump centre, while at larger radii, the
Stamatellos methods produces column density estimates that are too high.
Combining the two methods seems to produce estimates that are reasonably close
to that expected. In this case, the Modified Lombardi method produces
estimates very similar to those from the Combined method.

.

Figure 13: Column density of one of the clumps in the fragmenting high-mass
disc shown in the right-hand panel of Figure 11. The dashed line is the SPH
density integrated from the clump centre to $r_{\rm clump}=10$ AU. The symbols
are, again, a randomly selected sample of the column density estimates using
the Stamatellos method (yellow open circles), Lombardi method (red open
squares), Combined method (blue stars), and the Modified Lombardi method (open
green triangles). The Combined and Modified Lombardi methods give very similar
values near the centre of the clump and they agree with the Stamatellos
method. Figure 14: Cut through the disc along a radial line that includes the
clump shown in Figure 13 (the $M_{\rm disc}=0.25$ M⊙) showing the column
densities from the SPH simulation (black crosses) and the estimate from the
Stamatellos method (yellow open circles), the Lombardi method (red open
squares), the Combined method (blue stars), and the Modified Lombardi method
(open green triangles). For $r>50$ au, the points for the Modified Lombardi
method mostly lie on top of those of the Combined method.

Finally, Figure 14 shows a radial cut through the $M_{\rm disc}=0.25$ M⊙ disc
(Fig. 11 right panel) that includes the clump at $r\sim 65$ AU, the column
density of which was shown in Fig. 13. All the methods capture the clump but,
again, the Lombardi method over-estimates the column density near the clump
centre, while the Stamatellos method slightly over-estimates the column
density in the wings of the clump. The Combined and Modified Lombardi methods
both closely recover the peak of the clump and the clump profile.

##  5 Summary

We have described and tested two new methods for approximating radiative
cooling in hydrodynamical models in discs. The Stamatellos method, estimating
the optical depth via the gravitational potential, has been shown to
overestimate the column density, and therefore underestimate the cooling rate
in discs (Mercer et al., 2018). The Lombardi method, in which the column
density derives from the local quantity of the pressure gradient, provides a
more accurate estimate of the radiative cooling rate (Lombardi et al., 2015)
but we show that the column density estimates become unreliable close to the
disc midplane. These issues are resolved by combining the two methods, to
achieve improved estimates in regions of both high and low optical depth. This
method is particularly suited to spherical geometries. For disc geometries, we
present the Modified Lombardi method, in which the Lombardi estimate is
adjusted close to the disc midplane using an estimate of the scale height in a
self-gravitating disc. Whereas Mercer et al. (2018) show that the Lombardi and
Stamatellos methods both fail to provide accurate estimates of the optical
depth in a disc with clumps, the Modified Lombardi method closely recovers the
column density throughout a clumpy disc and traces the peaks and troughs
associated with clumps and spirals well. Of course, since the optical depth is
estimated from local disc properties, more complex geometries will not be
modelled well. Full radiative transfer is required to include the effects of
shadowing, for example.

It is straightforward to include the effects of a central star or stars,
external stellar companions and/or a variable interstellar radiation field
since the background temperature can incorporate arbitrary external heating.
Radiative transfer can be incorporated via the flux limited diffusion
approximation as formulated by Forgan et al. (2009), since the implementation
of the new methods is identical to that of Stamatellos et al. (2007). A future
development will be to consider the heating and gravitational influence of
multiple sink particles.

The Modified Lombardi method is ideal for modelling self-gravitating discs and
presents exciting prospects for studying young protoplanetary discs with a
physical treatment of cooling.

## Acknowledgements

We thank the anonymous reviewer for their kind comments that helped improve
the manuscript. This research made use of the DiRAC Data Intensive service at
Leicester, operated by the University of Leicester IT Services, which forms
part of the Science and Technology Facilities Council (STFC) DiRAC HPC
Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding
via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC
Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.
AKY and KR are grateful for support from the UK STFC via grant ST/V000594/1.
MC was supported by the Gates Cambridge Scholarship. RAB thanks the Royal
Society for their support via a University Research Fellowship. EC
acknowledges support from STFC grant ST/X508329/1 and DS from STFC grant
ST/Y002741/1. This work made use of numpy (Harris et al., 2020), matplotlib
(Hunter, J. D., 2007) and splash (Price, 2007). Fig. 4 made use of
WebPlotDigitizer (Rohatgi, 2022).

## Data Availability

The code developed for this paper and used to perform the simulations is
available from https://github.com/alisonkyoung1/phantom, which was forked from
phantom (https://github.com/danieljprice/phantom). The simulation setup files
are available on request.

## References

  * Ayliffe & Bate (2011) Ayliffe B. A., Bate M. R., 2011, MNRAS, 415, 576
  * Bell & Lin (1994) Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
  * Bertin & Lodato (1999) Bertin G., Lodato G., 1999, A&A, 350, 694
  * Booth & Clarke (2019) Booth R. A., Clarke C. J., 2019, MNRAS, 483, 3718
  * Boss (2017) Boss A. P., 2017, ApJ, 836, 53
  * Brandenburg & Das (2020) Brandenburg A., Das U., 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 162
  * Castor (2004) Castor J. I., 2004, Radiation Hydrodynamics
  * Clarke & Carswell (2014) Clarke C., Carswell B., 2014, Principles of Astrophysical Fluid Dynamics
  * Cossins et al. (2009) Cossins P., Lodato G., Clarke C. J., 2009, Monthly Notices of the Royal Astronomical Society, 393, 1157
  * Forgan et al. (2009) Forgan D., Rice K., Stamatellos D., Whitworth A., 2009, MNRAS, 394, 882
  * Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hunter, J. D. (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Johnson & Gammie (2003) Johnson B. M., Gammie C. F., 2003, ApJ, 597, 131
  * Laughlin & Bodenheimer (1994) Laughlin G., Bodenheimer P., 1994, ApJ, 436, 335
  * Levermore & Pomraning (1981) Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321
  * Lin & Kratter (2016) Lin M.-K., Kratter K. M., 2016, ApJ, 824, 91
  * Lodato & Rice (2004) Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630
  * Lombardi et al. (2015) Lombardi J. C., McInally W. G., Faber J. A., 2015, MNRAS, 447, 25
  * Mamatsashvili & Rice (2010) Mamatsashvili G. R., Rice W. K. M., 2010, MNRAS, 406, 2050
  * Masunaga & Inutsuka (2000) Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350
  * Mercer et al. (2018) Mercer A., Stamatellos D., Dunhill A., 2018, MNRAS, 478, 3478
  * Meru & Bate (2010) Meru F., Bate M. R., 2010, MNRAS, 406, 2279
  * Meru & Bate (2011) Meru F., Bate M. R., 2011, MNRAS, 411, L1
  * Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
  * Nealon et al. (2020) Nealon R., Price D. J., Pinte C., 2020, MNRAS, 493, L143
  * Nixon et al. (2018) Nixon C. J., King A. R., Pringle J. E., 2018, MNRAS, 477, 3273
  * Pinte et al. (2006) Pinte C., Ménard F., Duchêne G., Bastien P., 2006, A&A, 459, 797
  * Price (2007) Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159
  * Price et al. (2018) Price D. J., et al., 2018, Publications of the Astronomical Society of Australia, 35, e031
  * Rice et al. (2003) Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS, 339, 1025
  * Rohatgi (2022) Rohatgi A., 2022, Webplotdigitizer: Version 4.6, https://automeris.io/WebPlotDigitizer
  * Shi & Chiang (2014) Shi J.-M., Chiang E., 2014, ApJ, 789, 34
  * Stamatellos et al. (2007) Stamatellos D., Whitworth A. P., Bisbas T., Goodwin S., 2007, A&A, 475, 37
  * Vorobyov et al. (2020) Vorobyov E. I., Matsukoba R., Omukai, Kazuyuki Guedel, Manuel 2020, A&A, 638, A102
  * Wilkins & Clarke (2012) Wilkins D. R., Clarke C. J., 2012, MNRAS, 419, 3368


# Sound from rotors in non-uniform flow

Michael Carley
<EMAIL_ADDRESS>
###### Abstract

An analysis is presented for the evaluation of the acoustic field of a
rotating source in a non-uniform potential flow. Other than the restriction to
low flow Mach numbers, the method is exact and general. The variation in
radiation properties with source angle is handled by representation as a
Fourier series in source angle, giving rise to an asymmetrically varying
acoustic field evaluated by summation of the series. The method is used to
develop an exact solution for the model problem of a rotor operating near a
cylinder in cross-flow and sample calculations demonstrate the accuracy of the
technique when compared to full numerical evaluation. The calculations show
changes of greater than one decibel in the acoustic field due to flow non-
uniformity at a flow Mach number of 0.15, a typical speed for propeller
aircraft at takeoff.

##  1 Introduction

The prediction and analysis of sound from rotating sources is a long-standing
problem in acoustics, for its theoretical interest and for its importance in
many industrial applications. The problem can be viewed as made up of the
determination of the source terms, typically related to the fluid dynamics of
the system, and to the propagation of the resulting sound. In this paper, we
consider propagation of sound from a known source, on the assumption that the
source can be determined by some other means.

Even on linear theory, the calculation of the sound radiated by a rotating
source, here modeled as a sinusoidally-varying source distributed on a ring,
contains many features which complicate the calculation. The development of
prediction methods dates back some hundred years to the early studies of noise
from aircraft propellers [1]. In this work, the source is assumed to be at
rest and to radiate into a stationary fluid. The next extension to the theory
was the inclusion of flow effects in the form of a uniform flow parallel to
the axis of the source [2, 3], corresponding to a propeller in forward flight,
or to a propeller fixed in a mean flow, as in a wind tunnel. The introduction
of flow generates effects not present in the stationary fluid case, because of
changes in radiation efficiency caused by the increased speed of the source
relative to the fluid, and convective effects on retarded time. The acoustic
field, however, remains symmetric about the axis of rotation of the source,
because of the axial symmetry of the flow.

The next advance in the theoretical treatment of the problem was the
introduction of flow which is not symmetric about the source axis. For a
uniform flow his arises when a propeller is inclined at some angle to the
incoming flow, for example during take-off [4, 5, 6, 7]. In this case,
although the flow is uniform, because it is not aligned with the axis of
rotation of the source, the radiation efficiency varies as a function of
source angle. Early studies concentrated on the flow-generated unsteady source
terms, which were known to be efficient radiators of sound and only later were
asymmetric convection effects considered. Even if we consider only a steady
source, of constant strength in its frame of reference, this variation in
radiation efficiency gives rise to an asymmetric field. The asymmetry can be
modeled in a number of ways. In this paper, we adopt the viewpoint of Hanson
[6] who predicts the acoustic field as a Fourier series in azimuthal angle
about the source axis, with the source also decomposed into a Fourier series
or azimuthal modes. When there is no flow, or the flow is not inclined with
respect to the rotor axis, each source mode gives rise to one field mode of
the same order. When the flow is inclined with respect to the source axis, a
single source mode generates multiple field modes, whose magnitudes depend on
the source characteristics, and on the angle and speed of the flow.

In this paper, we extend the analysis of radiation from ring sources to the
case where the source radiates into a non-uniform flow, a generalization from
the existing theories for uniform flow. The only limitation is that the theory
used is correct to first order in flow Mach number. In underwater applications
this is not a significant restriction, but it does mean that the full range of
speeds relevant to aeronautical applications is not covered. In practice,
however, noise from aircraft is most troublesome during take-off and landing
when aircraft are near the ground, and in these cases aircraft are typically
operating at low flight speed, so the restriction to low Mach number is not as
limiting in practice as it might appear.

##  2 Analysis

Analysis of the problem of radiation from rotating, or ring, sources is
handled by treating the source as having sinusoidal variation in azimuthal
angle around the ring. Each such term gives rise to an integral which can be
expressed exactly by a series expansion and, if necessary, the contribution
from terms can be summed to give the overall field at any given frequency. In
this paper, we follow the same procedure, but we must take account of the
azimuthal variation in radiation properties caused by the non-uniform flow.
This is done by representing the variation in flow potential as a set of
azimuthal modes which convert the problem of radiation in non-uniform flow
into a sum of modes radiating into a static fluid.

In the rest of the paper, we do not consider the determination of the
magnitude of source terms, typically found from the aerodynamics of the
problem, nor do we extend the analysis to include interaction of the radiated
field with bodies in the flow. This is consistent with previous analyses [4,
5, 6] which have looked only at the effects of flow on direct radiation from
the source, without considering the secondary effect of scattering by nearby
bodies.

###  2.1 Sound from ring sources

Figure 1: Geometry and notation for radiation from ring source

The first step in evaluating the field from a ring source is a standard
integral which has been studied extensively for many years [2, 8, 4, 6, 9].
Figure 1 shows the geometry and notation. We adopt cylindrical coordinates
$(r,\psi,x)$ with the source lying on a ring of radius $a_{0}$ in the plane
$x=0$. The source has strength $s(\psi)\exp[-\mathrm{j}\omega t]$, with the
time variation suppressed from this point onwards. The acoustic potential at a
field point $(r,\psi,x)$ is given by integration on the source ring,

$\displaystyle\phi(r,\psi,x)$
$\displaystyle=\int_{0}^{2\pi}\frac{\mathrm{e}^{\mathrm{j}kR}}{4\pi
R}s(\psi_{1})\,\mathrm{d}\psi,$ $\displaystyle R^{2}$
$\displaystyle=r^{2}+x^{2}+a_{0}^{2}-2a_{0}r\cos(\psi-\psi_{1}),$

with $c$ speed of sound and wavenumber $k=\omega/c$. Subscript $1$ denotes a
variable of integration. Decomposing the source term into a Fourier series in
azimuthal angle $\psi_{1}$,

$\displaystyle s(\psi_{1})$
$\displaystyle=\sum_{n=-\infty}^{\infty}s_{n}\mathrm{e}^{\mathrm{j}n\psi_{1}},$

the potential is rewritten as a sum of terms

$\displaystyle\phi(r,\psi,x)$
$\displaystyle=\sum_{n=-\infty}^{\infty}s_{n}\mathrm{e}^{\mathrm{j}n\psi}I_{n}(k,a_{0};r,x).$

The integral $I_{n}(k,a_{0};r,x)$ is the acoustic potential radiated by a
single azimuthal mode of the source,

$\displaystyle I_{n}(k,a_{0};r,x)$
$\displaystyle=\int_{0}^{2\pi}\frac{\mathrm{e}^{\mathrm{j}(kR+n\psi_{1})}}{4\pi
R}\,\mathrm{d}\psi_{1},$ (1) $\displaystyle R^{2}$
$\displaystyle=x^{2}+r^{2}-2a_{0}r\cos\psi_{1}+a_{0}^{2},$

and is taken as the quantity of interest in the rest of the paper. The
integral has an exact series representation which will be used later to
calculate the acoustic field [9],

$\displaystyle I_{n}(k,a_{0};r,x)$
$\displaystyle=\frac{\mathrm{j}^{2n+1}}{2}\sum_{q=0}^{\infty}v_{n}h_{n+2q}^{(1)}(k\rho)j_{n+2q}(ka_{0})P_{n+2q}^{n}(x/\rho),$
(2) $\displaystyle v_{q}$
$\displaystyle=(-1)^{q}(2n+4q+1)\frac{(2q-1)!!}{(2n+2q)!!},$
$\displaystyle\rho^{2}$ $\displaystyle=x^{2}+r^{2},$

where $(\cdot)!!$ is the double factorial, $h_{n}^{(1)}(\cdot)$ is the
spherical Hankel function of the first kind, $j_{n}(\cdot)$ is the spherical
Bessel function, and $P_{n}^{m}(\cdot)$ is the associated Legendre function.

The form of the acoustic field generated by the ring source is controlled by
the parameter $M_{t}=ka_{0}/n$. For a steady source rotating in a uniform
flow, this is the Mach number associated with the rotation of the source, and
arises as a parameter in the analysis of propeller noise where it corresponds
to the velocity at the “tip” of a rotor blade. When asymmetry is introduced
into the system, this gives rise to modes of varying order $n$, at constant
wavenumber $k$. The parameter also arises in Tyler and Sofrin’s classic study
of propagation in ducts [10] where the effect of changing the mode order is
described in terms of changes to $M_{t}$. In particular, the propagation
properties change when $M_{t}>1$, the case of a “supersonic mode”. This
happens when a non-uniformity of some kind generates extra modes of order
$n\pm m$, $m=1,2,3,\ldots$. Should there be a value of $m$ such that $m\geq
n-ka_{0}$, $M_{t}\geq 1$ and the $n-m$ mode can radiate efficiently. If the
amplitude of this mode is large enough, this can lead to large changes in the
acoustic field.

The remainder of this paper extends the analysis for a ring source to the
problem of radiation into non-uniform flow by converting it into a sum of
integrals of the form of Equation 1, using two principal tools, Taylor’s
transformation for propagation in low Mach number potential flow, and
Panchekha’s formula for expansion of an exponential.

###  2.2 Sound propagation in non-uniform flows

Taylor’s transformation [11, 12, 13] is a method for converting problems of
propagation in non-uniform flows at low Mach number into equivalent problems
in a stationary fluid, by modifying the phase term of the Green’s function as
a function of the flow potential. The flow is specified as a velocity
potential $cM\Phi(\mathbf{x})$ where $M$ is the flow Mach number and
$\Phi(\mathbf{x})$ is the velocity potential for unit velocity. To first order
in Mach number, the acoustic potential at a point $\mathbf{x}$ due to a point
source at $\mathbf{x}_{0}$ is given by

$\displaystyle\phi(\mathbf{x})$
$\displaystyle=\frac{\mathrm{e}^{\mathrm{j}kR}}{4\pi
R}\mathrm{e}^{\mathrm{j}kM\Phi(\mathbf{x}_{0})}\mathrm{e}^{-\mathrm{j}kM\Phi(\mathbf{x})}.$
(3)

The prediction of sound in a low Mach number non-uniform flow thus becomes a
phase-shifted version of the problem in a stationary fluid. The difficulty of
the ring source problem lies in representing this phase shift in the radiation
integral of Equation 1 in a manner amenable to analysis.

###  2.3 Sound from ring sources in non-uniform flows

Inserting Taylor’s transformation, Equation 3, into Equation 1 yields an
integral representation for the acoustic potential in a general potential flow
at low Mach number. The integral can be evaluated numerically but this may be
time-consuming and yields no physical insight. We now show how the phase term
of the Taylor transformation can be represented as a set of azimuthal modes
converting the integral of Equation 1 into a sum of modes radiating into a
stationary fluid.

Upon substitution of the Taylor transformation into Equation 1, the radiation
integral becomes

$\displaystyle\int_{0}^{2\pi}\frac{\mathrm{e}^{\mathrm{j}(kR+n\psi_{1})}}{4\pi
R}\mathrm{e}^{\mathrm{j}kM\Phi(\psi_{1})}\,\mathrm{d}\psi_{1}.$ (4)

Our task now is to rewrite the exponential containing the flow potential in a
form suitable for analytical treatment.

The method used is based on the result of Panchekha [14] for series expansion
of the exponential of a Taylor series:

$\displaystyle\exp(a_{0}+a_{1}x+a_{2}x^{2}+\ldots)$
$\displaystyle=\mathrm{e}^{a_{0}}\sum_{n}A_{n}x^{n},$ (5) $\displaystyle
A_{n}$ $\displaystyle=\sum_{p\vdash n}\prod_{k,l\in p}\frac{a_{l}^{k}}{k!},$
(6)

where $p\vdash n$ denotes partitions $p$ of $n$, where $p$ is a set of terms
$k$ and $l$ such that $k_{1}l_{1}+k_{2}l_{2}+\ldots=n$.

In our case, we require an expansion for the exponential of a function
represented as a Fourier series. Writing $\mu=\exp[\mathrm{j}\psi]$,

$\displaystyle f(\psi)$
$\displaystyle=\sum_{q=-\infty}^{\infty}F_{q}\mathrm{e}^{\mathrm{j}q\psi}=\sum_{q=-\infty}^{\infty}F_{q}\mu^{q}.$
(7)

To develop a Fourier series which can be used for the Taylor transformation,
we write

$\displaystyle f(\psi)$ $\displaystyle=F_{0}+f_{+}(\mu)+f_{-}(\mu),$
$\displaystyle f_{+}(\mu)$ $\displaystyle=\sum_{q=1}^{\infty}F_{q}\mu^{q},$
$\displaystyle f_{-}(\mu)$ $\displaystyle=\sum_{q=1}^{\infty}F_{-q}\mu^{-q},$

and applying (5) yields

$\displaystyle\mathrm{e}^{\alpha f(\psi)}$ $\displaystyle=\mathrm{e}^{\alpha
F_{0}}\left(\sum_{q=0}^{\infty}A_{q}^{(+)}\mu^{q}\right)\left(\sum_{q=0}^{\infty}A_{q}^{(-)}\mu^{-q}\right),$
(8) $\displaystyle A_{q}^{(+)}$ $\displaystyle=\sum_{p\vdash q}\prod_{k,l\in
p}\frac{\left(\alpha F_{l}\right)^{k}}{k!},\quad A_{q}^{(-)}=\sum_{p\vdash
q}\prod_{k,l\in p}\frac{\left(\alpha F_{-l}\right)^{k}}{k!},$ (9)
$\displaystyle A_{0}^{(+)}$ $\displaystyle=A_{0}^{(-)}=0,$ (10)

and the Fourier expansion of the exponential is given by

$\displaystyle\mathrm{e}^{\alpha f(\psi)}$
$\displaystyle=\sum_{n=-\infty}^{\infty}c_{n}\mathrm{e}^{\mathrm{j}n\psi},$
(11) $\displaystyle c_{0}$ $\displaystyle=\mathrm{e}^{\alpha
F_{0}}\sum_{u=0}^{\infty}A_{u}^{(+)}A_{u}^{(-)},$ (12) $\displaystyle c_{q}$
$\displaystyle=\mathrm{e}^{\alpha
F_{0}}\sum_{u=1}^{\infty}A_{u+q}^{(+)}A_{u}^{(-)},$ (13) $\displaystyle
c_{-q}$ $\displaystyle=\mathrm{e}^{\alpha
F_{0}}\sum_{u=q+1}^{\infty}A_{u-q}^{(+)}A_{u}^{(-)}.$ (14)

Given a source of azimuthal order $n$, the acoustic potential at some field
point $\mathbf{x}$ is then

$\displaystyle\mathrm{e}^{-\mathrm{j}kM\Phi(\mathbf{x})}\mathrm{e}^{\mathrm{j}n\psi}\int_{0}^{2\pi}\frac{\mathrm{e}^{\mathrm{j}(kR+n\psi_{1})}}{4\pi
R}\mathrm{e}^{\mathrm{j}kM\Phi(\psi_{1})}\,\mathrm{d}\psi_{1}$
$\displaystyle=\mathrm{e}^{-\mathrm{j}kM\Phi(\mathbf{x})}\sum_{q=-\infty}^{\infty}c_{q}\mathrm{e}^{\mathrm{j}(n+q)\psi}\int_{0}^{2\pi}\frac{\mathrm{e}^{\mathrm{j}(kR+(n+q)\psi_{1})}}{4\pi
R}\,\mathrm{d}\psi_{1},$ (15)

where the coefficients $c_{n}$ are determined using the previous results with
$\alpha=\mathrm{j}kM$ and $F_{q}$ the coefficients of the Fourier expansion of
the flow potential on the source ring. The integrals on the right hand side
can be determined using the exact series expression of Equation 1.

###  2.4 Summary of method

The calculation method for the problem of a ring source of order $n$ radiating
into a non-uniform flow of Mach number $M$ can be summarized:

  1. 1.

compute the flow potential $\Phi(\mathbf{x})$ for unit velocity;

  2. 2.

compute the coefficients of the Fourier series of $\Phi(\mathbf{x})$ on the
source ring;

  3. 3.

apply the method of Section 2.3 to determine the coefficients $c_{q}$;

  4. 4.

use the series of Equation 2 and the summation of Equation 15 to determine the
acoustic field.

In the general case, the Fourier coefficients of $\Phi(\mathbf{x})$ may have
to be determined numerically, but in the next section we present a model
problem for which the required terms can be evaluated analytically.

##  3 Model problem

The method as presented so far is quite general as long as the Fourier
expansion of the flow potential can be determined. In practice, this will have
to be done numerically for many problems, but in the case of a source
interacting with a cylinder, a model problem for a rotor near a wing, the
Fourier coefficients can be found analytically. The problem contains the
essential features of the modification of an acoustic field by a realistic
potential flow and it will be seen that the flow gives rise to noticeable
effects on the overall field.

###  3.1 Geometry and notation

Figure 2: Geometry and notation

The model problem is shown in Figure 2. A wing is represented by an infinite
cylinder of radius $a$ lying along the $z$ axis. The source of radius $a_{0}$
is placed with its axis parallel to the $x$ axis and its center at
$(x_{0},y_{0})$. In calculating the flow potential, cylindrical coordinates
$(u,\theta,z)$ are used. In this coordinate system, the source center lies at
$(u_{0},\theta_{0})$ and points on the ring are given by
$(x_{0},y_{0}+a_{0}\sin\psi_{1},a_{0}\cos\psi_{1})$.

Figure 3: Streamlines for flow in positive $x$ direction around a cylinder of
radius $a$.

The solution for potential flow over a cylinder of radius $a$ can be found in
elementary fluid mechanics textbooks and in the notation of this paper is

$\displaystyle\Phi(\mathbf{x})$ $\displaystyle=x\left(1+a^{2}/u^{2}\right),$
(16)

for flow of unit velocity in the positive $x$ direction. Figure 3 shows the
streamlines for flow over a cylinder. The non-uniformity of the flow is clear
near the cylinder, and it is also clear that the flow is nearly uniform more
than about one cylinder radius from the cylinder surface.

In order to evaluate the effect of this flow on radiation from the ring
source, we expand $\Phi$ in a Fourier series in $\psi_{1}$, the angular
coordinate on the ring source of Figure 2:

$\displaystyle\Phi(\psi_{1})$
$\displaystyle=x_{0}\left(1+\frac{a^{2}}{u_{0}^{2}+2a_{0}y_{0}\sin\psi_{1}+a_{0}^{2}\sin^{2}\psi_{1}}\right),$
$\displaystyle=\sum_{q=-\infty}^{\infty}\Phi_{q}\mathrm{e}^{\mathrm{j}q\psi_{1}},$
$\displaystyle\Phi_{q}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\Phi(\psi_{1})\mathrm{e}^{-\mathrm{j}q\psi_{1}}\,\mathrm{d}\psi_{1}.$

Coefficients of the Fourier series can be evaluated using tabulated integrals
[15],

$\displaystyle\Phi_{0}$
$\displaystyle=x_{0}+\frac{a^{2}}{Ba_{0}}\cos(\gamma-\theta_{0}),$
$\displaystyle\Phi_{q}$
$\displaystyle=(-\mathrm{j})^{q}\frac{a^{2}}{Ba_{0}}\left(B^{2}-2BS\cos\gamma+S^{2}\right)^{q/2}\cos\left[(q-1)(\theta_{0}-\delta)+\gamma-\delta-q\pi/2\right],\,q>0,$
$\displaystyle\Phi_{-q}$ $\displaystyle=\Phi_{q}^{*},$

where

$\displaystyle B$ $\displaystyle=\left(S^{4}+2S^{2}\cos
2\theta_{0}+1\right)^{1/4},$ $\displaystyle S$ $\displaystyle=u_{0}/a_{0},$
$\displaystyle\gamma$ $\displaystyle=\frac{1}{2}\tan^{-1}\frac{\sin
2\theta_{0}}{S^{2}+2\cos 2\theta_{0}},$ $\displaystyle\delta$
$\displaystyle=\tan^{-1}\frac{B\sin\gamma}{B\cos\gamma-S}.$

The exponential containing the flow potential can then be written using
Equation 11 with $\alpha=\mathrm{j}kM$ and $F_{q}=\Phi_{q}$,

$\displaystyle\mathrm{e}^{\mathrm{j}kM\Phi(\psi_{1})}$
$\displaystyle=c_{0}+\sum_{q=1}^{\infty}\left(c_{q}\mathrm{e}^{\mathrm{j}q\psi_{1}}+c_{q}^{*}\mathrm{e}^{-\mathrm{j}q\psi_{1}}\right),$

with

$\displaystyle c_{q}$ $\displaystyle=\sum_{p=0}^{\infty}A_{p+q}A_{-p},$
$\displaystyle A_{n}$ $\displaystyle=\sum_{p\vdash n}\prod_{k,l\in
p}\frac{(\mathrm{j}kM\Phi_{l})^{k}}{k!},$ $\displaystyle A_{-n}$
$\displaystyle=\sum_{p\vdash n}\prod_{k,l\in
p}\frac{(\mathrm{j}kM\Phi^{*}_{l})^{k}}{k!}.$

The field can then be computed as

$\displaystyle\phi(\mathbf{x})$
$\displaystyle=\sum_{q=-\infty}^{\infty}c_{q}\mathrm{e}^{\mathrm{j}((n+q)\psi-
kM\Phi(\mathbf{x})}I_{n+q}(k,a_{0};r,x),$

reducing the general problem of radiation into a non-uniform flow to a
summation of elementary fields generated by individual azimuthal modes.

###  3.2 Results

To apply the results of the paper to a representative problem, parameters have
been estimated for the ATR class of turboprop aircraft. The propeller diameter
and speed are 3.93 m and 1200rpm respectively, so that the rotation Mach
number for the source is 0.73. The cylinder diameter has been set equal to
twice the maximum thickness of the ATR wing, with $a=0.29\mathrm{m}$. Flow
Mach number $M=0.15$, a typical value for takeoff of a turboprop aircraft.
Results have been calculated for a baseline case of $n=6$, representing the
lowest harmonic of the tonal noise from a six-bladed propeller. Figure 4 shows
the relative positions of the rotor and cylinder, with silhouettes of the
three configurations considered. The source is placed at negative $x_{0}$ to
correspond to a rotor ahead (upstream) of a wing leading edge.

|
---|---
Figure 4: Rotor positions for parametric study: left-hand side, notation;
right-hand side, relative positions of rotor and cylinder. Case a:
$(x_{0},y_{0})=(-a_{0},-a)$; b: $(x_{0},y_{0})=(-5a/2,-a_{0}/2)$; c:
$(x_{0},y_{0})=(-5a/4,-a_{0}/2)$.  Figure 5: Magnitude of computed acoustic
potential with $M=0.15$ flow, scaled on no-flow result. From top to bottom:
case a, b, c, and c with $n=8$. Solid line analytical result, dashed line
numerical evaluation.

Sample calculations for the problem are presented in Figure 5, in the form of
the magnitude of acoustic potential with flow, scaled on the no-flow case. As
a check on the method, the integral of Equation 4 was evaluated numerically
and its magnitude is also plotted. Results are shown for $n=6$ in each of the
configurations of Figure 4 and for $n=8$ in the third configuration. In each
case the field is calculated in the source plane on a circle of radius
$r=2a_{0}$. In considering the magnitude of flow effects on the radiated
noise, it is useful to note that Hanson found that a flow of Mach number
$M=0.2$ at 10 ° to the rotor axis altered the sound level by about 1.5 d B .
In our case, the flow is nominally at 0 ° and all asymmetry in the system
arises from the flow non-uniformity.

The first plot, for a rotor some distance upstream of the cylinder, shows a
very small change in the acoustic potential. This is as might be expected from
Figure 3 since the rotor lies in a region of flow which is almost uniform and
very little asymmetry will be introduced by the changes in potential around
the source ring. Moving the source closer to the cylinder, however, generates
quite large changes in the acoustic field which, we emphasize, are purely the
result of changes in the radiation properties of the source, which is held
constant throughout the calculations. In the second plot, with the source
brought within two radii of the cylinder surface, the field is modified by up
to 0.3 d B . In the third plot, moving the source closer still, a quarter
radius from the cylinder surface, generates a larger change in the acoustic
field, of some 0.6 d B . The largest change shown, however, comes with the
rotor in the same position, close to the cylinder, but with $n$ increased from
6 to 8, approximating the replacement of a six-bladed propeller with an eight-
bladed. The change in the acoustic field is now greater than 1 d B ,
comparable to the effect of a flow at an angle of 10 ° in Hanson’s study of
incidence effects. This effect is caused by the introduction of new modes with
higher $M_{t}$ radiating efficiently into the field.

##  4 Conclusions

An analysis has been developed for the prediction of sound from rotating
sources in a low Mach number non-uniform potential flow, giving an exact
formulation for the radiated field as a sum of azimuthal modes whose amplitude
is determined by the local flow. The theory has been applied to the model
problem of a rotor operating in the vicinity of a cylinder, for which the
required coefficients can be found analytically. Sample calculations compare
well to numerical evaluation of the acoustic field and the results show
changes in the field which vary strongly with rotor position relative to the
cylinder. The effect of flow non-uniformity on radiation properties has been
found to be comparable in magnitude to the effect of incidence in propeller
acoustics. Future work will consider disk sources, which have finite radial
extent modeling propeller blades of finite span, and the effects of wing lift
on the radiated noise, as finite circulation on the cylinder introduces
changes in the flow potential. An important open question remains that of the
relationship between the model problem and realistic systems which require
computational aeroacoustics methods for their prediction.

## References

  * [1] E. J. H. Lynam and H. A. Webb.  The emission of sound by airscrews.  Reports and memoranda 624, Aeronautical Research Committee, March 1919\.
  * [2] L. Gutin.  On the sound field of a rotating propeller.  Technical Memorandum 1195, NACA, Langley Aeronautical Laboratory, Langley Field, Va. USA, 1948.
  * [3] I. E. Garrick and C. E. Watkins.  A theoretical study of the effect of forward speed on the free-space sound-pressure field around propellers.  Report 1198, NACA, 1953.
  * [4] R. Stuff.  Noise field of a propeller with angular inflow.  AIAA Journal, 26(7):777–782, July 1988.
  * [5] R. Mani.  The radiation of sound from a propeller at angle of attack.  Proceedings of the Royal Society of London. A., 431:203–218, November 1990.
  * [6] D. B. Hanson.  Sound from a propeller at angle of attack: a new theoretical viewpoint.  Proceedings of the Royal Society of London. A., 449:315–328, 1995.
  * [7] L. M. B. C. Campos.  On the influence of angle-of-attack or sideslip on far-field propeller noise.  ASME Journal of Vibration and Acoustics, 121:50–58, 1999.
  * [8] S. E. Wright.  Sound radiation from a lifting rotor generated by asymmetric disk loading.  Journal of Sound and Vibration, 9(2):223–240, 1969.
  * [9] Michael Carley.  Series expansion for the sound field of a ring source.  Journal of the Acoustical Society of America, 128(6):3375–3380, 2010.
  * [10] J. M. Tyler and T. G. Sofrin.  Axial flow compressor noise studies.  Transactions of the Society of Automotive Engineers, 70:309–332, 1962.
  * [11] K. Taylor.  A transformation of the acoustic equation with implications for wind-tunnel and low-speed flight tests.  Proceedings of the Royal Society of London. A., 363:271–281, 1978.
  * [12] K. Taylor.  Acoustic generation by vibrating bodies in homentropic potential flow at low Mach number.  Journal of Sound and Vibration, 65(1):125–136, 1979.
  * [13] Anurag Agarwal and Ann P. Dowling.  Low-frequency acoustic shielding by the Silent Aircraft airframe.  AIAA Journal, 45(2):358–365, 2007.
  * [14] Pavel Panchekha.  Taylor expansions of Taylor expansions.  https://pavpanchekha.com/blog/taylor-exp.html, 2014.  Accessed 28 September 2017.
  * [15] I. Gradshteyn and I. M. Ryzhik.  Table of integrals, series, and products.  Academic, London, 5th edition, 1980.