relatively well but lacks the nonmonotonic behaviour of the simulated resonance frequency toward the high-power end of the spectrum. The CHRW approximation is nonmonotonic and is quite accurate when compared to the full numerical solution with $G/\omega_{\rm d} < 1$, and even for higher values of $G$ when $\omega_0/\omega_{\rm d} <1$ [@Lu12]. However, it still fails to match exactly the average numerical resonance at high drive powers, since it neglects the second and higher harmonics of the drive. At low powers, the deviation between the analytic and numerical resonance locations is caused by the vacuum BS shift $\chi_{\rm BS}=g^2/(\omega_0+\omega_{\rm c})=2\pi \times 0.7$ MHz. Strong driving clearly amplifies the BS effect, which appears as a deviation of nearly 3 MHz between the RWA and CHRW reflection minima.

From Fig. \[fig:BSshift\](b) we can see the existence of three regimes: at small average photon numbers, the response of the system is quantum-mechanical and corresponds to a cavity with constant frequency shift. When the number of photons becomes of the order of unity, the response is sensitive to the addition or removal of photons, indicating a dispersive photon blockade [@Houck2011]. At large number of photons, the photon blockade breaks [@Carmichael15], and the Bloch–Siegert shift 