# Integral of a harmonically phase-modulated sinusoid

I am trying to calculate integrals of the form $$I_n = \int_0^1 \cos\left[2\pi nt + \sum_{k=1}^N \beta_{nk} \sin(2\pi kt) \right] \mathrm{d}t,$$ where $n \in \{1, \ldots, N\}$ for some positive integer $N > 1$, and $\beta_{nk} \in \mathbb{R}$ for all $n,k \in \{1, \ldots, N\}$. In short, this question concerns whether there is an efficient way of calculating these integrals, with the help of analytics (series representation, etc.).

Such integrals arise, for example, when a set of oscillators, spread evenly in frequency, phase-modulate each other, and we aim to calculate the average position of each oscillator.

Clearly, these integrals can be computed independently for each $n$ by numerical quadrature, but I have found that the scaling is rather poor when the modulation indices $\beta_{nk}$ are not particularly small (say, $\beta_{nk} \sim 0.1$) and the number of oscillators $N$ is large (say, $N \sim 500$); presumably, this is caused by the rather high-frequency content of the integrand. I would like to know if anyone is aware of efficient ways to approximately compute the set of $I_n$—or even just $I_1$—either analytically (perhaps unlikely) or via an efficient series representation.

When the modulation indices are small, we can expand $I_n = I^{(0,1)}_n + \mathcal O(\beta^2)$, where the leading order term is \begin{align} I^{(0,1)}_n &= \int_0^1 \left[\cos(2\pi nt) - \sum_{k=1}^N \beta_{nk} \sin(2\pi nt)\sin(2\pi kt) \right] \mathrm{d}t \\ &= -\sum_{k=1}^N \beta_{nk} \int_0^1 \frac 1 2 \Big(\cos\big[2\pi(n-k)t\big] - \cos\big[2\pi(n+k)t\big]\Big); \end{align} that is, $I_n^{(0,1)} = -\beta_{nn}/2$. Of course, this approximation only holds when every $\beta_{nk}$ is small; otherwise, we need to take into account the terms at $\mathcal O(\beta^2)$.

I am also vaguely aware that phase-modulated integrals have a representation using Bessel functions. However, without having gone through all the details myself, it seems that at least the first-order term of the Bessel expansion produces a similar result to the above first-order Taylor expansion in $\beta$. Is there reason to believe that the Bessel expansion provides either more tractable higher-order terms, or that it converges faster? If so, what would be a good reference to consult in obtaining higher-order terms of a Bessel expansion?

Are there any other ways to represent the integral that would be more efficient than direct numerical quadrature, perhaps utilizing the somewhat special "up-to-$N$th harmonic" form of the phase modulation?