computing wavefront deviations for off-axis mirrors Suppose you have two off-axis aligned mirrors which are nearly planar, with small deformations such that a wavefront at plane $P_I$ at distance $-D$ from the mirrors, is refocused at distance $D$ at plane $P_O$. I only care about monochromatic waves of wavelength $\lambda$ 
I want to compute how much orders of Zernike mode expansion are required to account in the deformations of the mirrors, such that an input gaussian beam with intensity $I_I$ is projected to a gaussian beam with intensity $I_O$, such that $\frac{I_O}{I_I} > 1 - f$. That is, only a fraction $f$ of the beam is lost at most from the $T_{00}$ mode
I've been looking into the literature but most of it does the wavefront calculations for on-axis elements. I wanted to find some proper literature for the off-axis case that would suit the particulars of my problem
 A: The short answer (and likely one you're not going to like) to your question is that you're going to need as many Zernike terms as it is experimentally found are needed to model the aberration accurately. Zernike polynomials are normalized so that they contribute equally to the mean square phase error, and this latter, to the first order, is what sets the Strehl ratio, which is effectively what you're after here. The simplest analysis of the relation between aberration and phase is explored in sections 9.2 and 9.3 of Born and Wolf "Principles of Optics" (mine is the sixth edition); in 9.3 the Maréchal criterion is discussed, which is limit on the total mean square aberration that is Maréchal's definition of "diffraction limited performance". You may be able to narrow your description down to a few key Zernike polynomials, but there is no fundamental physics or mathematics to tell you which ones: either you must determine this experimentally or sometimes you can hazard a good guess if you can find any details of the production process that builds your mirrors out and understand the production process's "symmetries".
As for the ellipticity of your problem (wrought by the 45 degree beam folding in your "periscope" arrangement), I believe it can be handled by simple co-ordinate transformations as detailed below that make the mirror system nominally axisymmetric. 
So your Zernike analysis will be the "wonted" or "normal" one, but you transform your elliptical domain into a circular one with the $x$-direction shrink embodied in equations (2) and (3) below. Thus your situation will have the following special symmetry considerations:


*

*Spherical aberration terms become astigmatic for example, astigmatic terms will tend to become spherical - likewise, the symmetry class of all aberrations on the mirror surfaces will be changed by the shrink embodied by (2) and (3);

*Your mirrors will themselves are ellipsoidal rather than spherical (which I'm sure you understand) to yield the focusing power when 45 degree tilted rather than orthogonal to the optical axis. I'm pretty sure (without further analysis) the nominal contours on the ellipses are going to be of the form $\frac{1}{2} x_M^2 + y_M^2 = r^2$ (see (2) and (3) - so astigmatic Zernike terms are going to be especially important.

*There could be "azimuthal" misalignment between the two principal axes of the mirrors: i.e. the major and minor axes of one mirror may not be quite parallel to one another. This could lead to Zernike terms of high azimuthal symmetry classes: tetrafoil, heaxfoil and octofoil ($\cos$ or $\sin(4 \theta)$, $\cos$ or $\sin(6 \theta)$,  $\cos$ or $\sin(8 \theta)$;


Are you diamond turning these mirrors, btw? Diamond turning is extremely accurate, even in non axisymmetric components - the biggest errors are likely to be spherical (or spherical in the transformed co-ordinates in your case) so I'm guessing the major aberrations are going to arise from misalignements of the principal mirror axes.
Now we get onto how the analysis looks. Given $D \gg \lambda$, Fraunhofer diffraction applies, so the input field undergoes the equivalent of the following three steps: 


*

*Diffract to mirror system so that the wave has a diverging spherical wavefront of curvature radius $D$;

*Mirror system, through weak but intended ellipsoidal shape, imparts phase to the diffracting spherical wavefront so that now the wavefront has the conjugate phase (intended function, so that it becomes a converging rather than diverging wave) TOGETHER WITH an unintended phasing $\exp(i\,\phi(X,Y))$, where $\phi(X,Y)$ is the unintended aberration as a function of the rectangular mirror surface co-ordinates $X$ and $Y$ (to be defined so as to take account of the "folding" owing to the two 45 degree slanted mirrors);

*Diffraction "back" to the image plane. 


I am thus thinking of this problem as a more complicated version of a diffracting field bouncing off a spherical mirror aligned to the spherical wavefront so that the field focusses back at its initial positing in the ideal case. The transformation undergone by a field $\psi(x,y)$ in the object plane to reach the image plane is thus:
$$\psi \mapsto \mathfrak{F}^{-1}\, \exp(i\,\phi(X,Y))\, \mathfrak{F}\, \psi\quad\quad(1)$$
where $\mathfrak{F}$ is the unitary two dimensional Fourier transform. The inverse transform comes from the fact that my equivalent system imparts the nominal conjugate phase to the wavefront and sends it "backwards" to the image plane.
As long as we deal with measurements in the image plane alone, we cannot tell which of the pair any aberration comes from, so we may as well represent the combined effect of the "periscope" pair by the one aberration function $\phi(X,Y)$. This aberration is the thing whose effect you wish to quantify.
Now for the definitions of $X$ and $Y$, which is where the 45 degree tilting is accounted for:
$$X = \frac{k\,x_M}{\sqrt{2}\,D}\quad\quad(2)$$
$$Y = \frac{k\,y_M}{D}\quad\quad(3)$$
here $x_M$ and $y_M$ are the physical distances measured along the surface of the effective mirror, $i.e.$ imagine the mirrors flattened out taking away their by their nominal ellipsoidal curvature designed to offset the field curvature arising from the first diffraction. The residual deviations from this design goal are then added together to get the aberration function $\phi$. Here the $x$ direction is along the plane containing the system chief ray's nominal folded path: i.e. if the tube (line joining the two mirror centres) of the "periscope" is horizontal, then the $x$-axis is horizontal and likewise if the periscope tube is in any other direction. The square root of two factor accounts for the mirror tilt: the beam spreads out further on the mirror in the $x$ direction than it does in the $y$ owing to the effective mirror's intersecting the beam at the slant.
So now we are left with the effect of $\phi(X,Y)$. Let the input field be:
$$\psi_I\left(x,y\right) = \frac{1}{\sqrt{\pi}\,\sigma} \exp\left(-\frac{x^2+y^2}{2\,\sigma^2}\right)\quad\quad(4)$$
where $\sigma$ is the input field's spotsize. The first diffraction sends $\psi$ into the $X-Y$ space by:
$$\psi_I\left(x,y\right)\mapsto \Psi(X,Y)=\frac{1}{2\,\pi\,\sqrt{\pi}\,\sigma} \int\limits_{-\infty}^\infty \int\limits_{-\infty}^\infty \exp\left(-i \left(x\,X +y\,Y\right)\right)\,\exp\left(-\frac{x^2+y^2}{2\,\sigma^2}\right)\,\mathrm{d}x\,\mathrm{d}y = \frac{\sigma}{\sqrt{\pi}} \exp\left(-\frac{1}{2}\,\sigma^2\,\left(X^2+Y^2\right)\right)\quad\quad(5)$$
and now the output field is:
$$\psi_o\left(x,y\right)=\frac{\sigma}{2\,\pi\,\sqrt{\pi}} \int\limits_{-\infty}^\infty \int\limits_{-\infty}^\infty \exp\left(-i \left(x\,X +y\,Y\right)\right)\,\exp\left(-\frac{1}{2} \sigma^2 \left(X^2+Y^2\right)\right)\,\exp(i\,\phi(X,Y))\,\mathrm{d}X\,\mathrm{d}Y\quad\quad(6)$$
Now we can generalize the above result by including the effects of a defocus by noting that the diffraction through axial distance $\Delta z$ of a unidirectionally propagating scalar field fulfilling the Helmholtz equation undergoes the transformation defined by:
$$\psi \mapsto \mathfrak{F}^{-1} \,\exp\left(i\,\Delta z \sqrt{k^2 - X^2 - Y^2}\right)\,\mathfrak{F} \,\psi \approx e^{i\,k\,\Delta z}\, \mathfrak{F}^{-1}\, \exp\left(-i\,\Delta z \frac{X^2+Y^2}{2\,k}\right) \,\mathfrak{F}\, \psi\quad\quad(7)$$
where $\psi$ is the field at the input plane and is transformed to the field on a parallel plane a distance $\Delta z$ in the direction of the field's propagation. The first Fourier transform $\mathfrak{F}$ splits the field into plane wave components with $x$-wavenumber $X$ and $y$-wavenumber $Y$, then the phasing term sandwiched between the two Fourier transforms imparts the right phase delay for each plane wave component (if a plane wave has $x$ and $y$ wavevector components $X$ and $Y$, then the $z$ component of the wavevector must be $\sqrt{k^2 - X^2 - Y^2}$ where $k$ is the wavenumber), then the last inverse Fourier transform assembles all these delayed plane wave components into the diffracted field. The approximation assumes a low numerical aperture field, so that $X$ and $Y$ are small compared to the wavenumber $k$ and the plane wave components all make small angles with the axial direction.
So now we can combine (1) and (7) to get a version of (6) generalized to where the light field is transformed by the whole system, followed by a defocus of $\Delta z$:
$$\psi_o\left(x,y\right)=\frac{\sigma}{2\,\pi\,\sqrt{\pi}} \int\limits_{-\infty}^\infty \int\limits_{-\infty}^\infty \exp\left(i \left(\phi(X,Y) - x\,X -y\,Y - \Delta z \frac{X^2+Y^2}{2\,k}\right)\right)\,\exp\left(-\frac{1}{2} \sigma^2 \left(X^2+Y^2\right)\right)\,\mathrm{d}X\,\mathrm{d}Y\quad\quad(8)$$
The Strehl ratio is the intensity ratio of the above quantity to its aberration free value (which is $\frac{1}{\sqrt{\pi}\,\sigma}$), assuming the output field's peak amplitude is at position $(x, y)$. So now we look carefully at the total aberration term that lowers the peak amplitude:
$$\phi(X,Y) - x\,X -y\,Y - \Delta z \frac{X^2+Y^2}{2\,k}\quad\quad(9)$$
and rewrite that in polar co-ordinates $(\rho, \theta)$. To begin applying Zernikes to this aberration, we must define the "outer radius" in the normalized mirror co-ordinates $(X, Y)$; the Gaussian $\exp\left(-\frac{1}{2} \sigma^2 \left(X^2+Y^2\right)\right)$ is well "contained" within a radius of, say, $R = 4 / \sigma$; you can rework the following for any value of $R$, although I am suggesting a working value is going to be very like what I have just suggested. With polar co-ordinates normalized so that $X^2+Y^2 = R^2$ corresponds to the whole aperture and therefore to $\rho = 1$ (i.e. $\rho^2 R^2 = X^2+Y^2$), our total aberration is:
$$\phi(\rho, \theta) - x\,R\,\rho\,\cos\theta -y\,R\, \rho\,\sin\theta - \frac{\Delta z\, R^2}{4\,k} (2\rho^2 - 1) + \frac{\Delta z\, R^2}{4\,k}\quad\quad(10)$$
Here you see the two tilt and defocus Zernike functions $\rho\,\cos\theta$, $\rho\,\sin\theta$ and $2\rho^2 - 1$ and they correspond to sideways system misalignments (corresponding to tilts on your mirror, which I assume you have the freedom to impart for the sake of alignment) and small errors in the axial position of the output plane: again, I assume you have the freedom to impart small axial displacements to the mirror system in your alignment controls. The net aberration, after all these adjustments have been optimized, will be your mirror aberration $\phi(\rho, \theta)$ with the defocus and tilt Zernikes removed. Let $\tilde{\phi}(\rho, \theta)$ be the net aberration after tilt and defocus Zernikes have been removed. Then the Strehl ratio is:
$$\begin{array}{lcl}\mathcal{S} &=& \left(\frac{\sigma^2\,R^2}{2\,\pi} \int\limits_0^\infty\int\limits_0^{2\pi} \rho\,\exp\left(i \tilde{\phi}(\rho, \theta)\right)\,\exp\left(-\frac{R^2}{2} \sigma^2 \rho^2\right)\,\mathrm{d}\theta\,\mathrm{d}\rho\right)^2 \\
&\approx & \left(1 - \frac{\sigma^2\,R^2}{2\,\pi} \int\limits_0^R\int\limits_0^{2\pi} \rho\, \frac{\left(\tilde{\phi}(\rho, \theta)\right)^2}{2} \,\exp\left(-\frac{R^2}{2} \sigma^2 \rho^2\right)\,\mathrm{d}\theta\,\mathrm{d}\rho\right)^2\\
&\approx& 1 - \frac{\sigma^2\,R^2}{2\,\pi} \int\limits_0^R\int\limits_0^{2\pi} \left(\tilde{\phi}(\rho, \theta)\right)^2 \,\rho\,\exp\left(-\frac{R^2}{2} \sigma^2 \rho^2\right)\,\mathrm{d}\theta\,\mathrm{d}\rho\end{array}\quad\quad(10)$$
The last expression is the one I believe you are looking for to describe the degradation of the mode peak. The last two steps came from assuming small aberration so that we expand $\exp(i\tilde{\phi}(\rho,\theta))$ as a Taylor series to quadratic terms, and then to note that the integral of the linear Taylor term vanishes because we have stripped away the mean value of $\tilde{\phi}(\rho,\theta)$ when we removed the tilt and defocus Zernike terms. Recall that the constant $R$ is the domain radius in your normalized co-ordinates defined by (2) and (3) and so has the dimensions of inverse length.
This analysis is not too different from that in Section 9.3 of Born and Wolf. The main difference here is the co-ordinate transformations (2) and (3) which map the physical mirror surface co-ordinates into the normalized co-ordinates suitable for Zernike analysis, and the "lopsided" imparting of the $\sqrt{2}$ factor in (2) accounts for the tilted beam path in the "periscope", so you will need to transform your datasets by (2) and (3) before doing your Zernike analysis. Mathematica encodes Zernike radial functions by ZernikeR[n, m, r] (radial class $n$, azimuthal symmetry $m$, $r$ the independent variable). Most wavefront analysis softwares that come with interferometers can do the kind of linear transformation on datasets you want - the one I use is 4D Technology "4Sight". The software "Durango" has a free evaluation version but I haven't used it. Failing all that, I have my own C++ library for doing Zernike analysis so if you need to use this, contact me by my email on my user page. 
Lastly, for the sake of theoretical completeness, one strictly shouldn't use Zernike functions for situations where beams are Gaussian apodised like yours. One should use functions that are orthogonal with respect to the weight function $\rho \,\exp\left(-\frac{R^2}{2} \sigma^2 \rho^2\right)$; the Zernikes are orthogonal with respect to the weight function $\rho$. But this should not make a great deal of difference as long as you choose your "cutoff" radius $R$ well.
