  We present self-consistent triaxial stellar systems that have analytic
distribution functions (DFs) expressed in terms of the actions. These provide
triaxial density profiles with cores or cusps at the centre. They are the first
self-consistent triaxial models with analytic DFs suitable for modelling giant
ellipticals and dark haloes. Specifically, we study triaxial models that
reproduce the Hernquist profile from Williams & Evans (2015), as well as
flattened isochrones of the form proposed by Binney (2014). We explore the
kinematics and orbital structure of these models in some detail. The models
typically become more radially anisotropic on moving outwards, have velocity
ellipsoids aligned in Cartesian coordinates in the centre and aligned in
spherical polar coordinates in the outer parts.
  In projection, the ellipticity of the isophotes and the position angle of the
major axis of our models generally changes with radius. So, a natural
application is to elliptical galaxies that exhibit isophote twisting. As
triaxial St\"ackel models do not show isophote twists, our DFs are the first to
generate mass density distributions that do exhibit this phenomenon, typically
with a gradient of $\approx 10^\circ$/effective radius, which is comparable to
the data.
  Triaxiality is a natural consequence of models that are susceptible to the
radial orbit instability. We show how a family of spherical models with
anisotropy profiles that transition from isotropic at the centre to radially
anisotropic becomes unstable when the outer anisotropy is made sufficiently
radial. Models with a larger outer anisotropy can be constructed but are found
to be triaxial. We argue that the onset of the radial orbit instability can be
identified with the transition point when adiabatic relaxation yields strongly
triaxial rather than weakly spherical endpoints.
