  The squeezed-limit bispectrum, which is generated by nonlinear gravitational
evolution as well as inflationary physics, measures the correlation of three
wavenumbers, in the configuration where one wavenumber is much smaller than the
other two. Since the squeezed-limit bispectrum encodes the impact of a
large-scale fluctuation on the small-scale power spectrum, it can be understood
as how the small-scale power spectrum "responds" to the large-scale
fluctuation. Viewed in this way, the squeezed-limit bispectrum can be
calculated using the response approach even in the cases which do not submit to
perturbative treatment. To illustrate this point, we apply this approach to the
cross-correlation between the large-scale quasar density field and small-scale
Lyman-$\alpha$ forest flux power spectrum. In particular, using separate
universe simulations which implement changes in the large-scale density,
velocity gradient, and primordial power spectrum amplitude, we measure how the
Lyman-$\alpha$ forest flux power spectrum responds to the local,
long-wavelength quasar overdensity, and equivalently their squeezed-limit
bispectrum. We perform a Fisher forecast for the ability of future experiments
to constrain local non-Gaussianity using the bispectrum of quasars and the
Lyman-$\alpha$ forest. Combining with quasar and Lyman-$\alpha$ forest power
spectra to constrain the biases, we find that for DESI the expected $1-\sigma$
constraint is ${\rm err}[f_{\rm NL}]\sim60$. Ability for DESI to measure
$f_{\rm NL}$ through this channel is limited primarily by the aliasing and
instrumental noise of the Lyman-$\alpha$ forest flux power spectrum. The
combination of response approach and separate universe simulations provides a
novel technique to explore the constraints from the squeezed-limit bispectrum
between different observables.
