  Gravitational lensing is the deflection of light rays due to the gravity of
intervening masses. This phenomenon is observed in a variety of scales and
configurations, involving any non-uniform mass such as planets, stars,
galaxies, clusters of galaxies, and even the large scale structure of the
universe. Strong lensing occurs when the distortions are significant and
multiple images of the background source are observed. The lens objects must
align on the sky of order ~1 arcsecond for galaxy-galaxy lensing, or 10's of
arcseonds for cluster-galaxy lensing. As the discovery of lens systems has
grown to the low thousands, these systems have become pivotal for precision
measurements and addressing critical questions in astrophysics. Notably, they
facilitate the measurement of the Universe's expansion rate, dark matter,
supernovae, quasars, and the first stars among other topics. With future
surveys expected to discover hundreds of thousands of lensing systems, the
modelling and simulation of such systems must occur at orders of magnitude
larger scale then ever before. Here we present `caustics`, a Python package
designed to handle the extensive computational demands of modeling such a vast
number of lensing systems.
