A revisit of the tonal noise of small rotors

In this study, asymptotic analysis of the frequency-domain formulation to compute the tonal noise of the small rotors in the now ubiquitously multi-rotor powered drones is conducted. Simple scaling laws are proposed to evaluate the impacts of the influential parameters such as blade number, flow speed, rotation speed, unsteady motion, thrust and observer angle on the tonal noise. The rate of noise increment with thrust (or rotational speed) is determined by orders of blade passing frequency harmonics and the unsteady motion. The axial mean flow influence can be approximated by quadratic functions. At given thrust, the sound decreases rapidly with the radius and blade number as the surface pressure becomes less intensive. The higher tonal harmonics are significantly increased if unsteady motions, although of small-amplitude, are existed, as indicated by the defined sensitivity function, emphasizing that the unsteady motions should be avoided for quiet rotor designs. The scaling laws are examined by comparing with the full computations of the rotor noise using the frequency-domain method, the implementation of which has been validated by comparing with experiments. Good data collapse is obtained when the proposed scaling laws, which highlights the dominant influence of the design parameters, are incorporated.

SZ, 0000-0001-8235-0706 In this study, asymptotic analysis of the frequencydomain formulation to compute the tonal noise of the small rotors in the now ubiquitously multi-rotor powered drones is conducted. Simple scaling laws are proposed to evaluate the impacts of the influential parameters such as blade number, flow speed, rotation speed, unsteady motion, thrust and observer angle on the tonal noise. The rate of noise increment with thrust (or rotational speed) is determined by orders of blade passing frequency harmonics and the unsteady motion. The axial mean flow influence can be approximated by quadratic functions. At given thrust, the sound decreases rapidly with the radius and blade number as the surface pressure becomes less intensive. The higher tonal harmonics are significantly increased if unsteady motions, although of small-amplitude, are existed, as indicated by the defined sensitivity function, emphasizing that the unsteady motions should be avoided for quiet rotor designs. The scaling laws are examined by comparing with the full computations of the rotor noise using the frequencydomain method, the implementation of which has been validated by comparing with experiments. Good data collapse is obtained when the proposed scaling laws, which highlights the dominant influence of the design parameters, are incorporated.

Introduction
The multi-rotor powered flying vehicles such as drones are now increasingly ubiquitous [1]. The associated noise pollution is becoming a public concern since the vehicles are often operated in populated urban. The vehicles' weight, which can influence the noise level, varies with the widening drone applications from leisure to industrial purposes. The rotor noise contains rich tonal and broadband contents due to the periodic rotatory motions of the blades and the surrounding complex turbulent flows. Both tonal and broadband noise contents are harmful for human health. For example, the human auditory evoked response can be affected by the broadband noise [2], while the tonal noise can contribute to the annoyance as human is sensitive to the pitch characteristics [3]. In this study, we focus on the rotor tonal noise.
The rotor noise has been studied for a long time. In the 1930s, Gutin [4] proposed a rotor noise model based on the dipole noise produced by concentrated aerodynamic force. The theory was extended by Deming [5] by modelling the dipoles on the whole rotor disc, and the prediction results were compared with experiments. An early study to account for the effect of blade thickness was also done by Deming [6], and a correction for nonlinear thickness effect was proposed by Arnoldi [7]. The study by Garrick and Watkins [8] considered the impact of retarded time due to the distributed forces (sources). The pioneering investigations reflected the fact that an essential element in studying the rotor noise is to model the noise sources generated aerodynamically. With the establishment of the acoustic analogy theory, which was initially developed by Lighthill in the free space [9] and was extended to problems with static and moving boundaries by Curle [10] and Ffowcs-Williams and Hawkings (FW-H) [11] respectively, the sources of aerodynamic noise can be rigorously modelled. As many other aeroacoustic problems, the rotor noise can be divided into thickness noise, loading noise and quadrupole noise [12]. Therefore, the goal of rotor noise computation is converted to estimating the aerodynamic variables on the blade surfaces. The quadrupole volume source is often acoustically inefficient and therefore omitted at low Mach numbers. With the increasing power of computers and the advances in numerical simulations, a common approach is to resolve the near-field flows numerically. Then sound extrapolation methods [11,13,14] can be used to compute the far-field acoustics. However, high-fidelity rotor simulations are still too costly for parametric study, and (semi-) theoretical analysis is still needed to capture the main acoustic characteristics.
Another issue for the rotor noise study is due to the rotation motion, which is the major cause of tonal noise that presents as high peaks at the harmonics of the blade passing frequency (BPF) in the noise spectra. One approach for studying the rotation motion effect was to investigate the acoustic characteristics with moving unsteady flows and thereafter noise sources [15][16][17][18][19], in which the change of retarded time and the Doppler effect were considered. In an alternative approach, the temporal variation of the sound source due to the rotor motion was replaced by distributed sources with azimuthal phase variation. The rate of spatial phase variation has a significant impact on both the sound radiation efficiency and the directivity pattern [20]. Studies were conducted on the structure of the sound field based on the integral with respect to the phase angle [21]. Series expansion of the integration for efficient sound field computation with generic source distribution was conducted by Carley [22]. For rotors, the influence of fluctuated unsteady motions shall be highlighted. A representative unsteady motion of rotors happens when the rotor axis has a non-zero angle with the mean flow. The non-axially symmetric and cross-flow can cause unsteady flow, modulate the source strength, and impact on the sound radiation [23,24]. In the theoretical analysis by Wright [25], it was suggested that periodic force distributions with small amplitudes fluctuations could lead to significant noise generation. For small rotors, McKay and Kingan [26] measured the unsteady angular velocities and showed that the small velocity variation can significantly affect the noise at high-order harmonics of the BPF.
The difference in dealing with the rotation motion leads to two main approaches for practical rotor tonal noise prediction based on the time and frequency domains respectively. The timedomain method is based on Farassat Formulation 1A [13], which is an integral solution of the FW-H equation [11]. Usually, the time sequences of the near-field flow variables are needed, and the retarted time is estimated in the integral solutions. The method is often combined with numerical simulations and has been employed for rotor noise prediction [27]  the distribution of the aerodynamic noise produced on the blade surface is discontinuous, but repetitive, in the azimuthal direction at a periodic of 2π/B, where B is the blade number. The source locations are also varying periodically in time. Consequently, the equivalent sources at a different frequency (the BPF) can then be determined by the Fourier series expansion coefficients. The frequency-domain formulation for rotor noise was developed by Hanson [28][29][30] based on the generalized acoustic analogy of Goldstein [31] that considers the (uniform) background mean flow. The noise source computation based on local velocities and aerodynamic loadings was explicitly provided, e.g. by the mean surface modelling approach [30]. The time-domain method is able to capture the transient variation of the acoustic processes. By contrast, for the tonal noise caused by the periodic motion of the rotors, the flows as well as the noise sources are quasisteady such that the frequency-domain formulation can be more efficient as the time sequences of the input flow variables are not needed.
In the aeroacoustic studies, it is often of interest to understand the simple dependence of the noise radiation on the key influential parameters. With the widening drone applications, there is no doubt stringent noise regulations will be proposed. A reasonable regulation, likes those for other flying vehicles, shall be based on the trend of rotor noise emission with the weight (thrust), rotation speed and dimension. For example, for civil aircraft, the permitted maximum noise levels are varying with the maximum takeoff mass. Simple scaling laws of the sound pressure on the rotating speed and thrust can, therefore, provide insight into the noise characteristics, benefit the selection of design parameters and guide the drone noise certification. However, it is surprising that few investigations on this aspect were conducted for the rotor noise although theories were well developed. In the early rotor noise studies, Deming [5] measured the influence of the thrust and blade on the noise. Hubbard [32] measured the influence of the tip Mach number of the rotor on transport aeroplanes. Aravamudan et al. [33,34] proposed a M 6 t law of the tonal noise for the geometrically similar rotors, where M t is the rotation Mach number. In that study, the observer angles for the comparisons were different to ensure identical arguments (which also depend on M t ) in the Bessel functions. The previous investigations on large propellers and helicopter rotors often involved various working conditions and wide parameter ranges. For the small rotors (the radius of which is less than 0.2 m in this study) equipped on the drones, the properties of low Mach number and small dimension may ensure the accuracy in the approximations that can lead to simpler and more complete scaling laws.
In this work, scaling laws of the rotor tonal noise on the influential parameters such as rotation speed (or thrust), blade number and rotor radius will be developed based on the asymptotic analysis of the frequency-domain formula for rotor noise [30]. The influence of the unsteady motions on the noise signature is highlighted. The analyses are compared with the full computation, which has been validated by anechoic chamber experimental tests. The original contribution of the work lies on the development of the unsteady mode-dependent scaling laws for the tonal noise of small rotors. The following part of this article is organized as follows. Section 2 revisits the frequency-domain method, based on which the asymptotic analysis is conducted in §3 to propose the rotor noise scaling laws. Section 4 presents the summary.

Revisiting the frequency-domain formulation for rotor noise
In this section, the frequency-domain formulation of rotor noise proposed by Hanson and Parzych [30] is revisited. Simplified formula and major steps in the computation are presented, and the notations are introduced for the convenience of scaling law analysis in §3. The derivations start from the equations for the thickness and loading noise of a rotor and where u ∞ is the mean flow velocity in the axial direction (y 1 ), B is the blade number, m = μB is the sound at the μth BPF, ω = mΩ is the acoustic angular frequency and Ω is the rotation angular velocity. The variable α is the the initial phase of a source panel. G is the Green's function, and the argument is computed with respect to the source location (r, y 1 , φ + α) represented in the cylindrical coordinate. The observer in the spherical coordinate has the radius R, polar angle θ and azimuthal angle ϕ. In the Cartesian system, the locations of the observer and source are (x 1 , x 2 , x 3 ) = (R cos θ , R sin θ cos ϕ, R sin θ sin ϕ) and (y 1 , y 2 , y 3 ) = (y 1 , r cos(φ + α), r sin(φ + α)), respectively. A nomenclature of the notations in this study can be found in the supplementary material. Then, the Green's function is computed as G(x, y) = e ikσ /4π R, where k = ω/c 0 is the acoustic wave number, and the amplitude radius and phase radius are expressed as follows: where β 2 = 1 − M 2 and M = u ∞ /c 0 is the mean flow Mach number. The source term for the thickness noise is Q = v · n, which is the flux on the normal direction on the blade surface. The integral kernel of the loading noise is short for the computation that where F 1 , F r and F φ are the aerodynamic forces in the axial, radial and azimuthal directions, respectively. For the rotor blade, the projection of the normal vector n in the radial direction is often small. Therefore, it can be assumed that F r , which is the projection of the loading on the radial direction, is negligible. Considering that the flow contains unsteady but periodic components in the φ-direction, the source terms are represented as follows: where N 1 and N 2 are the number of the unsteady modes for velocity and loadings, respectively. The values are usually case dependent in practice. Then, equations (2.1) and (2.2) are reduced to and

(a) Modelling the source terms
The next step is to compute the source terms using the (known) aerodynamic variables [30]. Considering that the velocities contain different unsteady modes the source of thickness noise can be computed as follows: Q, u and v are summations of different orders of unsteady modes (described by exp(inφ) with n ranging from −N 1 to N 1 ), each term of the thickness noise source is computed as follows: For the rotor, the axial and azimuthal loading can be computed as where Φ is the sum of the advanced angle and the induced angle (see the supplemental material). F L = C L ρ 0 (u 2 + v 2 )/2 and F D = C D ρ 0 (u 2 + v 2 )/2 are the lift and drag, and C L and C D are the lift and drag coefficients. Recalling that the velocity contain different unsteady modes, the quadratic terms are then computed as follows: Under this assumption, it is known that N 2 = 2N 1 , and F n is determined as follows: Therefore, it is easy to know that the coefficient due to aerodynamic forces are In practice, in the presence of unsteady modes, the lift coefficient C L and drag coefficient C D might vary due to the unsteady aerodynamics effect [35]. Nevertheless, for simplicity, we assume that they are unchanged for the small-amplitude unsteady motions studied in this work.

(b) Comparison with experiments
To apply the frequency-domain formulation for noise prediction, the needed input flow variables can be computed using the blade element and moment theory (BEMT). In the blade element analysis, the lift and drag of each airfoil section are determined by the local velocities and the lift and drag coefficients C L and C D . In practice, C L and C D vary with the effective angle of attack (AoA), which is impacted by the local pitch angle and velocity components, and the local Reynolds number Re that is dependent on the local chord and velocity. In practical computations, the values are estimated based on a preconstructed C L and C D database of for airfoils at various AoA ∈ (−10 • , 20 • ) and Re ∈ (5 × 10 3 , 3 × 10 6 ) by a panel method [36]. The equations are enclosed by employing the moment theory that links the total momentum changes (represented by the velocity) with the lift and drag. The resultant differential and integral equations are solved iteratively to obtain the flow variables. The noise computations for several small-scale rotors are compared with the experimental tests on a rotor noise test rig placed in an anechoic chamber at the Hong Kong University of Science and Technology. The aerodynamic force is measured using a load cell, and the rotor noise is measured using a linear microphone array (with 10 GRAS 46BE free-field microphones) at 1.5 m from the rotor centre. The resultant observer angle ranges from θ ∈ (56 • , 118 • ), where θ = 0 • is the downstream direction. The schematic of the experiment is shown in figure 1a. The measured thrust of an APC 9x5 rotor at various revolutions per second (RPS), which were also measured by Brandt and Selig [37], are shown in figure 1b. The prediction of the thrust using BEMT (denoted by the curve) are close to both independent experiments (denoted by symbols). Then, the velocities, experiments (with good repeatability as indicated by the multiple, closely located symbols at each condition) at various observer angles and RPS, suggesting that the noise prediction model and its implementation are reliable. The experimental acoustic results shown in figure 2 are based on the noise data measured before the recirculation [38].

(c) Far-field approximation
Under the assumption that R r and R y 1 , the amplitude and phase radius can be computed are introduced for simplicity. Then, the Green's function can be approximated as follows: Then, it is straightforward to know ∂G/∂y 1 = −ikη 1 G and ∂G/∂φ = ikrη 2 sin(φ + α − ϕ)G. The thickness noise and loading noise can now be evaluated as follows: , of which the integral with respect to ψ ∈ (0, 2π ) is reduced to 2π (−i) m+n J m+n (kη 2 r) due to Bessel's first integral (cf. [39]) . Considering that sin ψ = (e iψ − e −iψ )/2i, we have and In practice, the remaining integral needs to be evaluated only once such that the computation time can be significantly reduced by the far-field approximation method.

Scaling laws of the rotor tonal noise
In practical applications, the rotor noise can be computed by using the frequency-domain formulation in §2, combined with a BEMT code. Systematic computation can be conducted to investigate the impact of the blade geometry for low-noise rotor design. However, it is still of interest to derive the simple dependence of the noise radiation on the key parameters such as flow speed, rotation speed, thrust and rotor dimension, despite that one may summarize the relations based on a number of full computation. Also, it is desired to highlight the impacts of different unsteady modes on the rotor tonal noise generation explicitly.
In this section, we analyse the rotor noise characteristics. Considering that for practical rotors in hovering and climbing, the amplitudes of the unsteady modes are often much smaller than (typically within 1%) the steady components [26]. Therefore, we can assume that the unsteady motions have small amplitudes compared with the steady components, i.e. |u n | |u 0 |, |v n | |v 0 |, for n = 0.
Then, high-order terms with the unsteady motion can be omitted in the analysis. Particularly, the term F n defined in equation (2.8) can be simplified as follows: where ζ n = 1 for n = 0, and ζ n = 2 otherwise. In addition, as shown in equation (2.12)  at μth BPF (such that m = μB) is determined by the Bessel function J m+n (krη 2 ). The meaning of different orders of unsteady modes is linked to the coefficients of the Fourier expansion of the non-uniform distribution of the sources in the φ (azimuthal) direction. For example, if the velocity varies sinusoidally as A cos(nφ) = A (exp(inφ) + exp(−inφ))/2, then the amplitudes for both the nth and −nth modes are A /2. In addition, considering that the Bessel function decays dramatically with n, at least for the small argument cases. Therefore, in the following discussions, we only consider the impact of the n ≤ 0 modes on the noise generation.
(a) Scaling laws of rotor noise in hovering Firstly, we study the rotors noise properties in hovering such that u ∞ = u 0 = 0. We denote the thickness noise, axial loading noise and azimuthal loading noise contributed by each modes n ≤ 0 as p T (x, n), p L1 (x, n) and p Lφ (x, n), respectively. For the small rotors (e.g. we focus on rotors with r < 0. where c (·) and r (·) are the integrals along the chord and radial directions respectively. Q n is influenced by both u n and v n , the amplitudes of which are usually unknown. However, we can assume they are proportional to the steady component v 0 = Ωr with case-dependent coefficients.
To propose a scaling law, we can represent the integral results by using tip radius r t . Considering that the chord of each section also increases with the tip radius, we can regard that the area is proportional to r 2 t . When performing the integral with respect to r, the kernel contains an exponent of r |m+n|+1 t , among which r |m+n| t is caused by J m+n (kη 2 r) and r t is due to the radial variation of Q n . Then, a coefficient of 1/(|m + n| + 2) will present in the results such that Furthermore, we denote κ n = (ĝ velocity distribution), the coefficients of 1/(|m + n| + 3) and 1/(|m + n| + 2) shall be included in the averaged quantity based on the tip radius r t . Then we have Ωr t 2c 0 |m+n|+3 sin |m+n| θ cos θ. (3.5) and Ωr t 2c 0 |m+n|+2 sin |m+n| θ, (3.6) where λ n and γ n , which have considered the sign due to ζ , denote the relative importance of the axial and rotational modes with the mean rotational speed. Based on equations (3.4), (3.5) and (3.6), we can analyse the scaling law properties of the rotor noise, the impact of different modes and directivity patterns, i.e. the dependences of the sound pressure on the observer angle θ . For the thickness noise p T and azimuthal loading noise p Lφ , the directivity patterns are determined by D 1 (θ ) = sin |m+n| θ. In case that |m + n| = 0, it has a maximum value at θ = π/2 and minimal values at θ = 0 (downstream) and π (upstream). For the noise contributed by the axial loading p L1 , the directivity function is D 2 (θ) = sin |m+n| θ cos θ. It is easy to see that D 2 (0) = D 2 (π/2) = D 2 (π ) = 0, suggesting that the corresponding noise has small amplitudes on the rotational plane and at the two poles. The derivative is D 2 (θ) = |m + n| sin |m+n|−1 θ cos 2 θ − sin |m+n|+1 θ . When m + n = 0 and sin θ = 0, there are maximum values of p L 1 at θ = cos −1 (±/ √ |m + n| + 1), where D 2 (θ ) = 0. However, for the unsteady modes that n = −m, the directivity function D 1 (θ) ≡ 1, suggesting that the noise emitted to different angles is identical. Considering that other modes have small quantities around θ = 0 and π , then the contribution by the n = −m mode can significantly increase the SPL at these angles. Also, for the axial loading p L1 , if n = −m, the directivity function is D 2 (θ ) = cos θ, which also has the maximum absolute quantities at θ = 0 and π . Figure 3 shows the dependence of the directivity functions with the unsteady mode n when m = 2. Also, as will be shown in the following text, the unsteady modes are acoustically efficient. Then, the modes with small amplitudes can lead to significant noise generation, causing high SPL distributions in the downstream direction that are often observed in experiments. . Therefore, the overall sound pressure may have dependence on M t with the power between |m + n| + 2 and |m + n| + 3, which is dependent on the order of the unsteady mode. Then, the sound pressure with the three different components can be evaluated as follows: where κ * n = κ n (|m + n| + 2)/(|m + n| + 3) and γ * n = γ n (m + n)/m. On the rotation plane, it is clear to see that p(θ = π/2) ∝ M |m+n|+2 t since θ = π/2. However, there is not a single scaling power for M t at all observer angles due to the presence of M t cos θ. To quantify the overall acoustic energy, we can define an averaged sound pressure as follows: For simplicity, we denote h 1 = κ * n + γ * n and h 2 = λ n such that the cos θ-dependent terms in p 2 can be written as h 2 1 + 2h 1 h 2 M t cos θ + h 2 2 M 2 t cos 2 θ. When performing integral with respect to θ ∈ (0, π ), the term α 2 M t cos θ vanishes because of cos θ is anti-symmetric with respect to θ = π/2. In addition, for the small rotors, we are studying, we assume M t 1 such that the h 2 2 M 2 t cos 2 θ can be omitted. Therefore, the averaged sound pressure due to the unsteady mode n could be scaled by M |m+n|+2 t , such that (3.9) where K n has included the impacts by the geometry and amplitudes of the unsteady motions and I ν can be evaluated using the recursive relation I ν = π 0 sin ν θ dθ = ν − 1 ν I ν−2 , with I 0 = π , and I 1 = 2. (3.10) It is also easy to know that I ν ≤ π for all ν. Equation (3.9) shows a simple dependence of p e on M t . Recalling that the thrust of a rotor can be computed as follows: where c T is thrust coefficient of a blade. The dependence of p e on T can be written as follows: where F(m, n) = I 1/2 2|m+n|+1 (μm/c T ) (|m+n|+2)/2 |m + n|!(|m + n| + 2)(2c 0 ) |m+n| . (3.12) Equation (3.12) also shows the dependence of p e on the unsteady mode amplitude K n , the thrust T, the rotor radius r t and the modes m and n. At a given thrust, p e decreases with r t since the aerodynamic loadings are less intensive for large areas.
In this work, we can examine the properties of the scaling laws proposed in this section by comparing with the full computation results. First, we focus on noise contribution by the steady mode of n = 0. In this case, p ∼ (mM t ) m+2 /(m!(m + 2)2 m ), which decays rapidly with the m = μB if mM t < 1. Therefore, the amplitudes for the higher-order harmonics of the BPF are very small such that we only focus on the tonal noise at BPF, i.e., μ = 1 and m = B. Figure 4 shows the results of different noise components for a two-bladed rotor with the rotation speed ranges from 60 to 120  Due to the influence of the axial loading noise, the total sound pressure is asymmetric about θ = π/2. The scaled results by Ω 4 also show relative good collapse.
The derivation also suggests that the sound pressure at the BPF (i.e. m = B) contributed by the n = 0 mode is significantly decreased with the blade number B due to the function F(m). Usually, for a given rotor geometry profile, c T can be approximately regarded as constant if the rotor is well designed. In this work, we compute the thrust of different rotor blades at various Ω. The rotors are geometrically similar, i.e. each rotor blade has the same twist angle and chord distribution c/r and the tip radius r t ranges from 0.08 to 0.17 m. The variation of the thrust T/B with Ω 2 r 4 is shown in figure 5a, where an approximately linear dependence is obtained. We also compute F(m, n) with m = B (at the BPF) and n = 0 (the steady mode), and the result is shown in figure 5b. The speed of sound is c 0 = 340 m s −1 , and c T = 1/50 is assumed. F(B, 0) decays with B rapidly, suggesting that the increment of blade can significantly reduce the rotor noise. In this work, we compute the averaged sound pressure contributed by the steady mode n = 0 of two-scaled rotors with r t ranges from 0.08 to 0.17 m and the rotation speeds ranges from 50 to 140 RPS. Consequently, the tip Mach number M t ranges from 0.074 to 0.44. The sound at the BPF of different blade numbers B = 2, 3 and 4 are studied. The sound pressure at the observer points with the radius of R = 15m are computed using the frequency-domain formulation [30] to estimate the averaged sound pressure p e defined in equation (3.9). The dependences of p e on F(B)T B/2+1 r −(B+1) t for the two scaled rotors are shown in figure 6, where F(B) = F(B, 0) is employed for simplicity. For each (scaled) rotor, the data are well collapsed as a linear function, suggesting that the dependence of the sound pressure on the thrust T and the tip radius r t is well captured. Also, for different B, the curves are close, which means that the influence of blade number on rotor noise is well estimated by F(B, 0). We have also tested various rotors with different geometries, and similar results were obtained. The slopes in the lines of the collapsed data indicate the acoustic performances of different rotors, i.e. the larger is the slope, and the more noise will it be at the given thrust.

(b) Impact of the unsteady modes
As denoted in equation (3.9), the noise contribution by each unsteady mode n is controlled by K n , which measures the ratio of the unsteady motions and v 0 = Ωr. This part of influence can be case dependent, and it should be determined as input for accurate prediction. Also, the acoustic efficiency of each unsteady mode is determined, which can be represented by the remaining part in equation (3.9). It is impacted by μ (order of BPF harmonics), n ≤ 0 (the unsteady mode) and M t that is the tip Mach number. In this work, we define a sensitivity function f (μ, n, M t ) = (μB) |μB+n|+2 |μB + n|!(|μB + n| + 2)2 |μB+n| I (3. 13) In case that n = 0, the function f (μ, 0, M t ) decays rapidly with μ, such that we only consider the rotor noise at BPF. By contrast, for the unsteady n < 0 modes, the acoustic properties can be significantly altered. As mentioned earlier, the unsteady mode can adjust the radiation pattern via  , which decays slower than m + 2 (i.e. the power for the steady mode).
To measure the importance of the unsteady mode, we study the ratio f (μ, n, M t )/f (1, 0, M t ), which quantifies the contribution by different modes to at different order of frequency compared to the usually studied noise at BPF mainly contributed by the steady aerodynamic forces. The results under different Mach numbers are shown in figure 7, in which f R (μ, n, M t ) = lg[f (μ, n, M t )/f (1, 0, M t )] is employed for better visualization. The larger the f R is, the easier the sound pressure can be produced at the corresponding BPF harmonic. For example, f R (μ = 3, n = −6, M t = 0.1) ≈ 3.2, which means that the n = −6 unsteady mode with its amplitude lower than Ωr/1000 (where 1000 is estimated from 10 f R ) can produce the sound at the third BPF harmonic, which has comparable SPL as the BPF contributed by the steady motion. The results for the n > 0 are also presented and small values are obtained, confirming that the noise produced by the positive modes is negligible. Approximately, it can be seen that the function f R reaches the peaks when n = −μB = −m. With the increase of M t , the function f R has large values in a wide range, as shown in figure 7b, which means that the unsteady modes can contribute to sound at different BPF harmonics more easily. The analysis confirms that the small amplitude of the unsteady loadings can lead to significant noise generation [25].
To highlight the influence of the unsteady motion, we compute the noise produced at different BPF harmonics (μ = 1, 2, . . .) contributed by different unsteady components. We assume that the amplitudes of 10 unsteady motions due to the fluctuation in rotation speed are v 0 /1000, which are sufficiently small. The observer distance is R = 2 m. The results at the observer angle θ = π/2 for two rotation speeds 60 and 120 RPS are shown in figure 8. It can be seen that each mode mainly contributes to a certain BPF harmonic, while the contributions on other harmonics decay rapidly. For both cases, it can be seen that the total sound contribution is dominated by the n = 0 mode (the solid blue line), suggesting that the use of steady aerodynamic forces can accurately capture at the BPF. However, for the tonal noise at higher BPF harmonics, the property can be significantly influenced by the although small-amplitude unsteady motions because of the rapidly increasing sensitivity function f (μ, n, M t ). As a result, the levels of the high-order BPF harmonics are comparable with the tonal noise at the BPF, which is supposed to be highest in most cases. The directivity patterns of the BPF harmonics m = μB contributed by the unsteady motions n = 1, . . . , 10 are shown in figure 9. It can be seen that the directivity patterns are significantly varied due to the difference in D 1 (θ ) = sin |m+n| θ (and D 2 (θ) = sin |m+n| θ cos θ) coupling with f (μ, n, M t ) at different m and n. For the sound at the BPF, the SPL at θ = 0 and 180 • , which has minimum values for cases of steady motion, is increased due to the unsteady motions. Also, the amplitude of the higher BPF components is increased and may have the maximum values in the upstream  and downstream directions. The analysis in this study confirmed that the unsteady motions are the main cause of the high-order BPF harmonics and can significantly alter the directivity, which is consistent with the analysis of the sensitivity function in equation (3.13).
To further show the effect of unsteady motions on the rotor noise generation, we consider a practical rotor with unsteady fluctuations. In practice, the small-amplitude unsteady rotations is difficult to measure experimentally, while the aerodynamic forces can be easily measured by using high-accuracy load cells (with a sampling frequency of 20 kHz). Therefore, we may approximately estimate the velocity fluctuation from the measured thrust signals by assuming T ∝ u 2 + v 2 . By definition, the velocities contain the unsteady components that u = u 0 + u n exp(inΩt) and v = v 0 + v n exp(inΩt). Considering that u 0 = 0, |u n | v 0 and |v n | v 0 , the quadratic terms of v n are omitted such that By performing a Fourier transform of the measured thrust T(t), we can obtain the amplitudes of the thrust at different frequencies nΩ. Then, we can approximately determine the amplitude of v n by v n /v 0 ≈ T n /2T 0 , where T n is the thrust amplitude at nΩ and T 0 is the mean thrust. Figure 10 shows the time-domain signal and the thrust spectrum of a rotor at 60 RPS measured in the experiment shown in figure 1. The acoustic and thrust data are acquired after the recirculation happens [38]. The spectrum is normalized by the zero-frequency T 0 (the mean thrust) to show the relative amplitude. It is seen that the amplitudes of the unsteady thrust are lower than 1% of the mean thrust. As analysed earlier, the relative values are used to measure the unsteady velocity as input for noise computation. The phases of the unsteady modes are configured as random. The predicted tonal noise at different BPF components is shown in . The difference of the observer distance R at different observer angle θ is due to the linear microphone array shown in figure 1. The results by the experimental measurements (red dashed lines), predictions with the unsteady modes estimated from the thrust (solid blue lines with symbols) and the prediction only considering the steady mode (solid green lines with symbols) are presented. At the BPF, both predictions with and without the unsteady modes are close to the experiment result. At higher orders of BPF, the prediction only considering the steady motion decrease rapidly. By contrast, with the incorporation of the unsteady modes, some of the highfrequency BPF components are overpredicted, especially at θ = 90 • . The possible reason is that the estimation of the velocity function is simply obtained by assuming that T ∝ v 2 , which may lead to the overestimation of the velocity fluctuation since unsteady aerodynamic effects are not considered [35]. The thrust fluctuation may also be caused by other effects that should not be attributed to the variation of velocity fluctuations. Nevertheless, the prediction informed by the measured unsteady thrust yields better agreement with the experiment. Also, it highlights that the unsteady motion leads to significant tonal noise production at high BPF harmonics. that approximations can be made by omitting the terms with orders higher than M 2 . Therefore, we have R 0 ≈ R(1 − (1/2)M 2 sin 2 θ ) such that Also, it is easy to see that η 1 ≈ cos θ + M + cos θ (1 + sin 2 θ/2)M 2 , so that We also have η 2 ≈ sin θ (1 + (sin 2 θ/2)M 2 ) such that krη 2 = mΩr sin θ/c 0 (1 + (sin 2 θ/2)M 2 ). Therefore, the Bessel functions are approximated as follows: (3.14) Being similar to the analysis for the hovering cases, for the thickness noise contributed by each mode n ≤ 0, the source strength is proportional to the combination of u n and v n . We denote the ratio κ n = Q n /v 0 to measure the relative importance of the unsteady fluctuations. Then, the thickness noise can estimated (based on the tip radius r t ) as follows: Ωr t 2c 0 |m+n|+2 sin |m+n| θg 1 (M, θ), (3.15) where g 1 (M, θ) = 1 + M cos θ + [1 + sin 2 θ/2(|m + n| + 1)]M 2 shows the modulation of Mach number on the directivity pattern. Similarly, for the loading noise, the source strengths due to the velocities are due to the quadratic terms. Then, we denote the ratios with (Ωr) 2 as λ n and γ n , respectively, such that the axial and azimuth loading noise components are evaluated as follows: Ωr t 2c 0 |m+n|+3 sin |m+n| θg 2 (M, θ) (3.16) and p Lφ (x, n) ≈ γ n ρ 0 c 2 0 πμR r t m |m+n|+1 (m + n) |m + n|!(|m + n| + 2) Ωr t 2c 0 |m+n|+2 sin |m+n| θg 3 (M, θ), (3.17) where g 2 and g 3 are evaluated as follows: . The axial mean flow can also alter the noise generation by impacting the coefficients κ n , λ n and γ n , which are characterized by the velocity amplitudes u n and v n . To further quantify the impact, we introduce the following notations v n = k n v 0 , and u n = α n k n u n , (3.20) where v 0 = Ωr is the rotation speed at the tip, and the coefficients α n to measure the relative amplitudes of the axial and azimuthal velocities are case dependent.  The unsteady modes n < 0 can first influence the sensitivity function f (μ, n, M t ) defined in equation (3.13). The impact is the same as for the hovering case that a small amplitude fluctuation can contribute to high SPL for higher BPF contents. The unsteady modes can also impact directivity patterns. The influence is on the order of power in sin |m+n| θ, which may lead to omnidirectional distributions when n = −m. They can also impact on the functions g 1 , g 2 and  Another important impact is that they can modulate the strength in κ n . In the climbing configuration, the streamwise fluctuations u n are coupled with the axial flow u 0 and the sound radiation is therefore dependent on J. The relative importance of each mode is also affected by the rotor geometry. For practical applications, the climbing speed of the rotors might be low such that G(J) may tend to a constant. However, it should be carefully noted that the effect of the even small aerodynamic mean flow cannot be regarded as negligible because unsteady aerodynamic motions can be likely produced. Then significant noise generation can appear due to the sensitivity function f (μ, n, M t ), which has tremendous values for the unsteady modes. The values of the unsteady aerodynamic motions can be case dependent. For quiet rotor designs, it is needed to avoid the occurrence of the unsteady motions [25].

Summary
This work studies the tonal noise characteristics of the small rotors. Based on the frequencydomain method of rotor noise [30], simple scaling laws of the noise dependence on rotation speed, thrust, unsteady modes, blade number, rotor size and axial flow are proposed. Mainly, the dependence of the sound pressure on the rotation Mach number M t is influenced by the order of BPF harmonic μ, the unsteady mode n and the blade number B as p ∝ M |μB+n|+2 t , i.e. the Mach number dependence is not governed by the conventional M 6 t law for aerodynamic noise. The results can also be represented by the thrust T as p ∝ T |μB+n|/2+1 , which can benefit the quick noise evaluation and to make reasonable noise regulations for drones equipped with small rotors. At a given thrust, the sound pressure depends on the blade radius by r −|μB+n|−1 t , suggesting that small radius will lead to intensive pressure distribution and therefore noise radiation. Also, the sound pressure decreases with the blade number rapidly as suggested by the function F(m = μB, 0). The impact of unsteady modes is quantified by the sensitivity function f (μ, n, M t ), which has a huge value for higher n modes. The sound properties can be significantly altered when unsteady periodic fluctuations are presented, especially at high BPF harmonics. A fluctuation of the rotational speed within 1% can lead to high BPF noise contents with comparable SPL as that at the BPF. A rotor with good acoustic design shall avoid the appearance of unsteady fluctuations. As for the impact of axial flow, which can be described by the advance ratio J, it can slightly modulate the directivity pattern. The influence on sound level can be characterized by a quadratic function of J, with the case-dependent coefficients determined by the rotor geometry. Also, we have conducted systematic full computations of the rotor noise, and the accuracy has been examined by experiments in anechoic chamber. The good collapse of the data are obtained when the proposed scaling laws are considered. Ethics. The manuscript is not currently being considered for publication in another journal. Data accessibility. The code and example input files can be found at http://aantc.ust.hk/Team/zhongsiyang/ sub_mainpage_tone.html.