  We use the IllustrisTNG suite of cosmological simulations to measure
intrinsic alignment (IA) bispectra of dark matter subhalos between redshifts 0
and 1. We decompose the intrinsic shear field into E- and B-modes and find that
the bispectra $B_{\delta\delta\mathrm{E}}$ and $B_{\delta\mathrm{EE}}$, between
the matter overdensity field, $\delta$, and the E-mode field, are detected with
high significance. We also model the IA bispectra analytically using a method
consistent with the two-point non-linear alignment model. We use this model and
the simulation measurements to infer the intrinsic alignment amplitude
$A_\mathrm{IA}$ and find that values of $A_\mathrm{IA}$ obtained from IA power
spectra and bispectra agree well at scales up to $k_\mathrm{max}= 2 \, h
\mathrm{Mpc}^{-1}$, for example at $z=1$ $A_\mathrm{IA} = 2.13 \pm$ 0.02 from
the cross power spectrum between the matter overdensity and E-mode fields and
$A_\mathrm{IA} =2.11 \pm$ 0.03 from $B_{\delta\delta \mathrm{E}}$. This
demonstrates that a single physically motivated model can jointly model
two-point and three-point statistics of intrinsic alignments, thus enabling a
cleaner separation between intrinsic alignments and cosmological weak lensing
signals.
