# Current status of PAPYRUS : the pyramid based adaptive optics system at
LAM/OHP

Muslimov E.c,Levraud N.a,b,Chambouleyron V.a,b ,
Boudjema I.a ,Lau A.a, Caillat A.a,Pedreros F.a,
Otten G.a,El Hadi K.a, Joaquina K.a,Lopez M.a,
El Morsy M.a,Beltramo Martin O.a,Fétick R.a,b,
Ke Z.a,Sauvage J-F.a,b, Neichel B.a, Fusco T.a,b,
Schmitt J.d, Le Van Suu A.d,Charton J.e,Schimpf A.e,
Martin B.e,Dintrono F.e, Esposito S.f,Pina E.f
a–Laboratoire d’Astrophysique de Marseille, France
b–DOTA,ONERA, Université Paris Saclay, F-91123 Palaiseau, France
c–NOVA Optical IR Instrumentation Group, Dwingeloo, Netherlands
d –Observatoire de Haute Provence, Saint-Michel-L’observatoire, France
e –ALPAO, Montbonnot-Saint-Martin, France
f –Osservatorio Astrofisico di Arcetri, INAF, Italie

###### Abstract

The Provence Adaptive optics Pyramid Run System (PAPYRUS) is a pyramid-based
Adaptive Optics (AO) system that will be installed at the Coude focus of the
1.52m telescope (T152) at the Observatoire de Haute Provence (OHP). The
project is being developed by PhD students and Postdocs across France with
support from staff members consolidating the existing expertise and hardware
into an R&D testbed. This testbed allows us to run various pyramid wavefront
sensing (WFS) control algorithms on-sky and experiment on new concepts for
wavefront control with additional benefit from the high number of available
nights at this telescope. It will also function as a teaching tool for
students during the planned AO summer school at OHP. To our knowledge, this is
one of the first pedagogic pyramid-based AO systems on-sky. The key components
of PAPYRUS are a $17\times 17$ actuators Alpao deformable mirror with a Alpao
RTC, a very low noise camera OCAM2k, and a 4-faces glass pyramid. PAPYRUS is
designed in order to be a simple and modular system to explore wavefront
control with a pyramid WFS on sky. We present an overview of PAPYRUS, a
description of the opto-mechanical design and the current status of the
project.

Keywords— Adaptive Optics, Pyramid sensor

##  1 INTRODUCTION

###  1.1 Presentation of the project

The Provence Adaptive-optics PYramid RUn System (PAPYRUS) is an adaptive
optics system to be installed at the Coude focus of the T152 telescope
(diameter of 1.52 m) at Observatoire de Haute Provence (OHP). The main
specificity of Papyrus is to use a pyramid wavefront sensor, which shows
better SNR performances than the usual Shack-Hartmann wavefront sensor. The
pyramid wavefront sensor has widely been described theoretically [1, 2] and
tested on bench [3]. PAPYRUS is a demonstrator of feasibility on sky and a TRL
maturation step in between the current tests performed at LAM and the future
generation of pyramid wavefront sensors such as ELT/HARMONI.

The specificity of PAPYRUS is to be a low cost and student led project. Indeed
PAPYRUS will use components already available at LAM and ONERA in order to
drastically reduce hardware costs. The PAPYRUS project encompasses all the
following steps:

  * •

System specification and requirements

  * •

System optical and mechanical design

  * •

Assembly and testing at LAM

  * •

Installation and testing at OHP

  * •

Maintenance and eventual upgrades at OHP

###  1.2 PAPYRUS project objectives

As stated above, PAPYRUS is a demonstrator of capability to place a pyramid
wavefront sensor on sky. The LAM/GRD and ONERA teams will gain a valuable
knowledge on the pyramid behaviour, especially regarding the development of
specific control laws. The impact of the PSF modulation on the pyramid will be
a major topic of study thanks to a dedicated modulation mirror and to a
dedicated pyramid focal plane camera [4] (“gain scheduling camera”).
Management of the pyramid optical gains and non-linearities is still a subject
under active research. PAPYRUS might be upgraded in the future to test new
techniques or components for astronomy with AO.

Moreover, PAPYRUS data will be used to check image post-processing
possibilities from pyramid AO systems. PSF analysis, PSF estimation [5] or
deconvolution [6] is envisaged.

Finally, PAPYRUS will be a pedagogical on-sky AO bench located near Marseille
(France) to teach on AO techniques students coming from different laboratories
or summer schools at OHP. Observing proposals may be issued by astronomers
once the system performances and sky coverage have been validated on sky.

###  1.3 The Observatoire de Haute Provence (OHP) site

The astronomical site of OHP, located inthe South-East of France, has been
chosen for accessibility reasons. It is the closest observatory to the
development laboratory at LAM, the OHP T152 telescope has available space for
installing an AO bench and available observation time. The telescope
observation time will have to be shared with the AURELIE instrument, the
latter being used only 50 (clear) nights per year approximately. The OHP
suffers from approximately 60-80 cloudy nights per year (figure 3), it remains
around 220-240 available clear nights per year for PAPYRUS.

The OHP is supporting the PAPYRUS effort. OHP will allow installation of the
AO bench on the telescope and teaching some PAPYRUS people to use the
telescope in autonomy. The T152 telescope is made of four mirrors of 88%
reflectivity each, providing a total transmission of 52%. Recent coating on
the primary mirror might slightly improve this transmission factor. The main
drawback is the T152 important pointing error of $\pm 2$ arcmin. Once the star
has been located in the T152 within this error range, a manual search would
have to be performed to make it appear on the PAPYRUS camera. This manual
searching phase will reduce slightly the effective AO lock time of PAPYRUS but
is perfectly manageable according to OHP astronomers and our own usage of the
telescope.

##  2 Adaptive Optics Design

###  2.1 Atmospheric conditions at OHP

The adaptive optics design and performances depends on the atmospheric
conditions at the observing site. Data from OHP [7] in Fig.1 show $\sim 60$
nights of good seeing $s<2"$. The median of the data is in between $s=2"$ and
$s=4"$, representing approximately $200$ nights per year. These values give a
typical Fried parameter $4<r_{0}<6$ cm at a wavelength of $600$ nm.

The median wind speed at OHP is about $5$ m/s, with bursts up to $7$ m/s.

Figure 1: Statistics of seeing at OHP in 2018. Image taken and adpted from OHP
website.

###  2.2 Adaptive optics error budget

The performances of the system may be stated in term of Strehl ratio between
the AO corrected long exposure and the T152 diffraction limit. Using the
Marechal approximation, this Strehl ratio is related to the residual
electromagnetic phase variance over the pupil as

$S_{r}=e^{-\sigma^{2}}$ (1)

The usual way of balancing an AO error budget[8] is to write the phase
variance $\sigma^{2}$ as the quadratic sum of different contributions

$\sigma^{2}=\sigma_{\text{fitting}}^{2}+\sigma_{\text{aliasing}}^{2}+\sigma_{\text{noise}}^{2}+\sigma_{\text{temporal}}^{2}+\sigma_{\text{others}}^{2}$
(2)

The definition and computations of these terms is given below. These
expressions being chromatic, we choose to give all values in the visible, at a
wavelength of $635$ nm.

####  2.2.1 Fitting error

The so-called fitting error is due to the limited number of actuators of the
deformable mirror. Thus the high spatial frequencies of the turbulent phase
cannot be corrected by the AO system and are present in the residuals. The
expression of this error is

$\sigma_{\text{fitting}}^{2}=0.27\left(\frac{d}{r_{0}}\right)^{5/3}$ (3)

with $d$ the pitch between two actuators and $r_{0}$ the Fried parameter of
the turbulence. For our $17\times 17$ deformable mirror on a $D=1.52$ m pupil,
the pitch is $9.5$ cm, giving an error of

$\sigma_{\text{fitting}}^{2}=0.6\text{\;rad}^{2}$ (4)

####  2.2.2 Aliasing error

The wavefront measurement is sampled for detection, producing aliasing of the
high frequencies onto the low frequencies. This information from the high
frequencies enters in the AO loop and generates a deformable mirror response.
However the aliasing is smaller on a pyramid wavefront sensor than ona Shack-
Hartmann wavefront sensor. It is then neglected in our AO budget

$\sigma_{\text{aliasing}}^{2}=0\text{\;rad}^{2}$ (5)

####  2.2.3 Noise error

The photon noise and detector noise both produce intensity fluctuations on the
wavefront sensor and consequently in the estimated wavefront. The advantage of
our AO design is to use a low-noise OCAM2K detector for wavefront sensing.
Detector noise is thus neglected.

Since PAPYRUS will be targeting bright stars in its early design, photon noise
will be mitigated. In order to achieve a correct SNR $=10$ on each
measurement, one needs $100$ photons per point of measure per frame. This flux
on the detector translates into flux of the AO guiding target as

$N_{\text{ph/m2/s}}=N_{\text{ph/meas/frame}}\frac{F}{\eta}\left(\frac{n_{meas}}{D}\right)^{2}$
(6)

where $N_{\text{ph/m2/s}}$ is the flux of the guiding target in photons per m2
per second, $N_{\text{ph/meas/frame}}=100$ is the desired flux on the detector
in photons per point of measurement per frame, $F=300$ Hz is the frame
frequency, $\eta\simeq 10\%$ is the global photon efficiency from telescope to
the WFS detector, $n_{meas}=16$ is the number of measurements in the diameter,
$D=1.52$ m is the telescope diameter. One finds

$N_{\text{ph/m2/s}}=33\times 10^{6}$ (7)

For a wavelength of $700$ nm and a filter width of $20$ nm, this flux is
obtained for stars of magnitude $V\leq 4$. With this condition on the star
magnitude the noise error propagated in the loop is neglected

$\sigma_{\text{noise}}^{2}=0\text{\;rad}^{2}$ (8)

####  2.2.4 Temporal error

The temporal error comes from the delay between the actual phase present on
the pupil and its correction by the deformable mirror. The temporal error
writes

$\sigma_{\text{temporal}}^{2}=0.04\left(\frac{V}{DB}\right)^{2}\left(\frac{D}{r_{0}}\right)^{5/3}\sum_{n=1}^{N_{r}}(n+1)^{-2/3}$
(9)

with $V$ the turbulence equivalent wind speed (Taylor frozen flow hypothesis),
$B$ the AO loop bandwidth and $N_{r}$ the number of radial modes corrected.
For a control law including an integrator, the equivalent bandwidth writes

$B=\frac{F}{2\pi}\sqrt{\frac{g}{1+2\tau F}}$ (10)

with $g$ the integrator gain and $\tau$ the loop delay. A reasonable choice of
gain ensuring both performances and robustness is

$g=\frac{1}{1+\tau F}$ (11)

Considering $\tau F=2$ frames delay for PAPYRUS and a wind speed $V=7$ m/s, we
deduce the bandwidth and then the temporal error

$\sigma_{\text{temporal}}^{2}=5.9\text{\;rad}^{2}$ (12)

####  2.2.5 Other sources of error

Other sources of error include, but are not limited to, scintillation on the
pupil, non-common path errors, and calibration errors. These errors are not
taken into account, especially in front of the important temporal error.

###  2.3 Conclusions on the budget

The error budget above is highly dominated by the temporal error. This error
is established with the highest wind speed at OHP, and thus corresponds to a
worst case of PAPYRUS usage. It is planned to mitigate this error by
implementing more complex control laws to increase the bandwidth, and using a
better RTC. Indeed we plan to replace our current computer (run on Matlab) by
a dedicated ALPAO high-speed RTC.

Adaptive-optics corrected PSF are simulated using OOMAO[9] end-to-end
simulations with the PAPYRUS system parameters, see Fig. 2. The temporal error
being the dominant one, we use three different AO loop frequencies for the
simulations. It shows that the current $F=300$ Hz is quite limitating and
incline us to go for a better RTC at $F=500$ Hz at least.

Figure 2: PSF from end-to-end simulations. Three cases of loop frequency are
considered to simulate the PAPYRUS PSF.

##  3 Optical Design

###  3.1 Initial data and design challenges

The optical design is built to work together with the T152 telescope at OHP.
It represents an F/28.56 Cassegrain system with the focal length of 43.4 m and
the maximum field of view (FoV) of $\pm 3^{\circ}$. The telescope has residual
aberrations, which don’t exceed $16.8\mu m$,while the Airy disk radius is
$13.95\mu m$.

Also, we had certain cameras available for the scientific and wavefront
sensing branches as well as the commercial deformable and tip-tilt mirrors.
Their parameters were also driving the design.

Thus, except of the telescope parameters the initial data for the optical
design development included:

  1. 1.

Diameter of off-the-shelf mirror used for slow tip/tilt =25.4 mm;

  2. 2.

Diameter of the existing DM = 37.5 mm;

  3. 3.

Max diameter of the existing fast tip/tilt mirror = 12 mm;

  4. 4.

Format of ORCA science camera = $13.312\times 13.312mm^{2}$ with $2048\times
2048$ pixels;

  5. 5.

Format of OCAM WFS camera= $5.76\times 5.7mm^{2}$ with $240\times 240$ pixels;

  6. 6.

The measured parameters of existing glass pyramid: facet angle =
$8.9^{\circ}$, material LF5, leading to deflection angle of $\pm
5.44^{\circ}$;

  7. 7.

Science camera FoV = $\pm 1"$;

  8. 8.

WFS FoV = $\pm 0.25:$;

  9. 9.

Working wavebands = 400-1000 nm (Science camera) and 630-690 nm (WFS with
filter).

In addition to these initial data there were some limitations, which made the
development challenging. First, the project budget and the development time
were limited. Therefore we had to rely only on off-the-shelf components and
iteratively repeat the design to find compromise solutions. Second, since our
demonstrator is targeting multiple tasks, its’ design should be flexible with
a possibility to replace some components, to connect testing and alignment
modules or to feed other instruments through optical ports, to get an access
to any par of it for educational purposes. Third, the optical design of the
bench should provide a high image quality to minimise the instrumental WFE and
the difference between the science and WFS branches. Finally, there were
limitations in overall dimension and optical interfaces because the bench will
be installed on a telescope among the existing instruments.

###  3.2 Optical design overview and analysis

With all of the above-listed initial values and limitations we developed the
following optical design (see Figure 3). In order to provide flexible
connection with the telescope and compensate the height difference the movable
periscope unit 2, 3 is used to pick up the beam from the telescope focal plane
1. Hereafter we use two optical relays to form intermediate focal and pupil
planes and place key components there. The relays are built on the basis of
off-axis parabolic (OAP) mirrors. use of OAP allows to exclude chrmatic
aberrations and significantly reduce the overall dimensions, although makes
the alignment more difficult. The beam is collimated (4) and the pupil is
imaged on the slow tip-tilt mirror 5. Further the beam is focused again (6)
with $-0.5^{x}$ magnification and the pupil is re-imaged to the infinity.
Folding mirror 7 represents an optical port, which can be used to couple the
telescope simulator optical to the intermediate focus 8. The beam is
collimated again (9) and the pupil is imaged onto DM 10. Then it is focused
(11) in the same way with $–1^{x}$ magnification and the pupil projection to
the infinity. The image formed in focal plane is detected by science camera
13. Before the science camera a beamsplitter 12 is mounted to feed the WFS
branch. Note, that the WFS branch uses lens optics, because the components
there should have short focal lenghts and work with relatively large angular
fields. The WFS $1^{st}$ intermediate focus 14 can be used to set a diaphragm.
The beam is re-collimated (15) and the pupil is projected to the modulation
mirror 17. In the formed parallel beam after 15 a filter 16 can be mounted to
minimize the WFS chromatism. It should be noted that some chromatic aberration
is inherent to the pyramid itself regardless of the auxiliary optics
properties. The beam is focused (18), while the pupil is projected to the
infinity. The focused beam is split (19) to form two identical images at the
tracking camera 20 and the tip of pyramid 21. The pyramid splits the beam in
quadrants with $5.44^{\circ}$ deflection. The resultant beams are collimated
(22) and the pupil is projected to the WFS camera 23.

Figure 3: General view of the PAPYRUS optical design (see the text for the
elements notation).

Thus, the optical design consists of a series of functional elements, each
served with a collimating and focusing relay. Since each focusing component
should project the pupil to the infinity, the pupils coincide with the
focusing components focal planes. For the same reason, the collimating
components form the pupil images in their back focal planes. The actual focal
lengths and magnifications in the system were choosen to fit the initial
design conditions and to obtain the focal lengths, $F/\\#$ and field values
close to those of available commercial components.

The Figure 4 shows the image quality in the field of view of science camera.
The spot diagram root mean square (RMS) radii vary between 5.3 and 14.9 $\mu
m$, while the Airy disk radius is $10.1\mu m$. So, the image quality is
diffraction limited on-axis and slightly degrades towards the field edge. The
WFE RMS value is $0.048\lambda$ with the maximum of $0.16\lambda$ at 700 nm.
The on sky quality will however be dominated by AO residuals and
anisoplanatism, so this design is completely accepted for PAPYRUS.

Figure 4: Image quality at the scientific camera: A – spot diagrams over the
FoV (the colors correspond to wavelengths in $\mu m$); B – wavefront error in
waves<EMAIL_ADDRESS>m$ referred to the exit pupil plane.

The image quality assessment for the WFS branch is shown in Figure 5. The
geometrical spot size remains stable across the field, since the nominal
pyramid FoV is only $\pm 14.4"$. It has the RMS radius of 5.1-6.0 $\mu m$ ,
while the Airy radius is 12.1 $\mu m$. The WFE for a single pyramid facet is
$0.046\lambda$ RMS at 700 nm. It becomes possible to reach the diffraction
limit only with use of a red band filter. Otherwise the chromatic aberrations
inherent to this system would blur the spot images by factor of 4.

Figure 5: WFS branch image quality at the tracking camera: A – spot diagrams
over the FoV (the colors correspond to wavelengths in $\mu m$ with filter); B
– wavefront error in waves<EMAIL_ADDRESS>m$ referred to the exit pupil plane.

We would like to note separately the design issue related to the last WFS
lens, which forms the image on the camera after the pyramid. It should work
with a relatively low $F/\\#$ and extended field of view, but also in an
unordinary setup, where the pupil is shifted and doesn’t coincide with the
nominal aperture stop or its’ images. In the meantime this lens should be an
off-the-shelf component and provide a high image quality. A $Canon^{TM}$ EF-S
24mm F/2.8 STM pancake photographic lens was chosen, since it is has small
total length and provide high image quality. However, to account for possible
vignetting and aberrations growth is this setup we had to reconstruct an
approximate optical design of the lens. An approximate algorithm of this
reverse-engineering is given below. We hope that it can be useful for those,
who are going to use commercial lenses in non-standard laboratory setups.

  1. 1.

The general data about the lens like the number of lenses and groups as well
as the composition of positive and negative components were taken from the
manufacturer[10].

  2. 2.

The closest patent was found[11]. The lens composition is similar to the known
one and the publication date is the closest ot the lens commercial release.

  3. 3.

The patent lens was scaled to fit the focal length and the first lens diameter
of the commercial lens.

  4. 4.

The lens was optimized to fit such criteria as the MTF values and vignetting
known from open sources. The radii of curvature and thicknesses were used as
free variables. In addition to standard boundary conditions defining the edge
and center thicknesses, the overall system length and diameters and optical
powers of individual components were restricted. The lens design after
optimization is shown in Figure 6,

Figure 6: Reconstructed design of the WFS lens: A – general view, B – MTF in
comparison with the reference curve from open sources.

The figure 7 shows the opto-mechanical implementation of PAPYRUS with the
different elements already presented before. The optomechanical design uses
off-the-shelf parts for all the adjustments and holding the components and has
as less as possible custom parts only as mechanical interface adapters. The
mechnaical mounts provide the necessary number of degrees of freedom for every
unit. The entire bench is assembled on a single breadboard. The optical path
is divided into straight segments and all the components of such a segment are
mounted on a single optical rail. Position of each of the optical rails is
given by a couple (line+corner) flat templates. The periscope assembly and the
cameras are mounted via custom adapters, which can be 3d-printed or milled in
aluminium. All the OAP mirrors have adapters to compensate the vertex lateral
shift. Finally, the WFS camera together with the pyramid are mounted in a
custom housing with adjustment rings and attached directly to the WFS camera
C-mount. This optomechanical implementation allows some room for future
upgades of the AO bench, and visiting scientific instrument to be adapted to
the bench.

Figure 7: Opto-mechanical model of PAPYRUS in 3D view (left) and top view
(right). GSC: gain scheduling camera. DM: deformable mirror. OAP: off-axis
parabola. Slow TT: slow tip-tilt mirror.

##  4 Current Status of papyrus

###  4.1 Assembly Integration and Tests

The AIT of PAPYRUS will be organized following this strategy:

  * •

component validation at LAM

  * •

global bench alignement and validation at LAM

  * •

final coupling and commissioning at T152 OHP.

####  4.1.1 Test of the pyramid

The critical components of PAPYRUS have been tested individually. The glass
pyramid is obviously one of these critical components. Pictures of the pyramid
tip (Fig. 8) show a ”rooftop” configuration due to manufacture precision. The
rooftop size is approximately the diffraction size that would prevent WFS
sensing in one direction if non modulated. However the modulation of the PSF
around the pyramid tip solves this issue.

Figure 8: Picture of the PAPYRUS pyramid tip, showing a ”rooftop”
configuration.

###  4.2 Testing the WFS branch

The WFS branch is the core of the AO instrument. This branch has already been
aligned at LAM (Fig. 9) to test the critical components to work together: the
pyramid, the WFS camera, the DM and the modulation mirror. We processed the
raw intensities and were able to compute the calibration matrix of the system
by sending push-pull commands on the DM.

Once the tests performed, the WFS branch has been unmounted. The alignment of
the full bench is currently under process at LAM.

Figure 9: Tests of the WFS branch at LAM. Left: picture of the WFS branch.
Middle: image seen on the WFS camera when a poke is sent to the DM. Right:
computed slopes maps.

###  4.3 On-sky measurements at OHP

Measurements at the Coude focus of the T152 telescope have been performed on
July 20th and 21st, 2021. The objective of these measurements is to
characterise the telescope pupil plane and focal plane. The possible high-
amplitude movements of these planes is problematic for the AO system. However
the pupil plane has shown stability better than $1\%$ of the diameter during a
long-exposure observation. The movement of the pupil is consequently less than
a fraction of point of measurement of the phase (equivalent to a sub-aperture
for a Shack-Hartmann WFS).

The star position is targeted within a $\pm 2$ arcmin uncertainty diameter in
the sky due to the telescope pointing lack of accuracy. However once the star
is found it can be correctly centred manually. Then the PSF suffers a small
drift of $\sim 20$ arcsec/hour. This drift is supposed to come from the
telescope weight unbalance that might accelerate or slow down the Earth
rotation compensation. However this drift is small and will be managed by the
PAPYRUS slow tip-tilt mirror.

Figure 10: Characterisation of the OHP/T152 telescope. Left: image of the
pupil. The dashed red ellipse shows presence of localised amplitude
fluctuations at the edge of the pupil. Right: long-exposure PSF obtained on
Vega.

###  4.4 Schedule

The schedule of the PAPYRUS project is given on Fig. 11. Integration at OHP is
planned for end 2021, and PAPYRUS will be available for the astronomical
community beginning or mid 2022.

Figure 11: Schedule of the PAPYRUS project.

##  5 Conclusion

As it was presented here PAPYRUS will soon allow LAM to improved its
operational expertise on Pyramids and deliver an AO corrected image at OHP.
The system has been designed, the components have been characterised, and the
first integration is underway with preparations for arrival at OHP also well
advanced. The project will soon be opened to the community and we will welcome
any new idea and collaboration to add new capabilities to the bench.

##  6 Acknowledgement

This project has been supported by ONERA, LAM and ALPAO in particular the WOLF
ANR. We also thanks First Light and the OHP for their contribution this
project.

## References

  * [1] O. Fauvarque, B. Neichel, T. Fusco, J.-F. Sauvage, and O. Girault, “Generalformalism for fourier-based wave front sensing: application to the pyramidwave front sensors,”Journal of Astronomical Telescopes, Instruments, andSystems3(1), p. 019001, 2017.
  * [2] O. Fauvarque, P. Janin-Potiron, C. Correia, Y. Brule, B. Neichel, V. Cham-bouleyron, J.-F. Sauvage, and T. Fusco, “Kernel formalism applied to fourier-based wave-front sensing in presence of residual phases,”JOSA A36(7),pp. 1241–1251, 2019.
  * [3] P. Janin-Potiron, V. Chambouleyron, L. Schatz, O. Fauvarque, C. Z. Bond,Y. Abautret, E. R. Muslimov, K. El Hadi, J.-F. Sauvage, K. Dohlen,et al.,“Adaptive optics with programmable fourier-based wavefront sensors: a spatial light modulator approach to the lam/onera on-sky pyramid sensor testbed,”Journal of Astronomical Telescopes, Instruments, and Systems5(3),p. 039001, 2019.
  * [4] V. Chambouleyron, O. Fauvarque, J.-F. Sauvage, B. Neichel, and T. Fusco,“The focal-plane assisted pyramid wavefront sensor: enabling frame-by-frame optical gains tracking,”arXiv preprint arXiv:2103.02297,2021.
  * [5] R. J. Fetick, L. Mugnier, T. Fusco, and B. Neichel, “Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach,”MonthlyNotices of the Royal Astronomical Society496(4), pp. 4209–4220, 2020.
  * [6] L. M. Mugnier, T. Fusco, and J.-M. Conan, “Mistral: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images,”JOSA A21(10), pp. 1841–1854, 2004.
  * [7] OHP, “meteorological data.” http://www.obs-hp.fr/climatologie.shtml.
  * [8] F. J. Rigaut, J.-P. Veran, and O. Lai, “Analytical model for Shack-Hartmann-based adaptive optics systems,” in Adaptive Optical System Technologies,3353, pp. 1038–1048, International Society for Optics and Photonics, 1998.
  * [9] R. Conan and C. Correia, “Object-oriented matlab adaptive optics toolbox,”inAdaptive optics systems IV,9148, p. 91486C, International Society for Optics and Photonics, 2014.
  * [10] Canon Inc.,“Canoncameramuseum.”https://global.canon/en/c-museum/product/ef434.html.
  * [11] M.Ohtake and M. Mori, “Compact wide-angle objective lens,” Sept. 16 1997.US Patent 5,668,669A


Lyapunov Function for the Nonlinear Moog Voltage Controlled Filter

Technical Report AAG2021a
Acoustics and Audio Group, University of Edinburgh

Stefan Bilbao

March 19 2021

In this short report, a new Lyapunov function for the Moog voltage-controlled
filter is demonstrated, under zero-input conditions, and under nonlinear
autonomous conditions (i.e. when parameters are not time-varying). The new
definition allows for a proof of stability over the entire allowable range of
parameters (cutoff frequency and resonance), and can be used as a starting
point for Hamiltonian-based numerical simulation methods.

##  1 Introduction: The Moog voltage-controlled filter

The operation of the Moog ladder filter or voltage-controlled filter (VCF)
[1], including nonlinear effects, is covered by Hélie [2, 3]. The digital
emulation of the Moog VCF dates back to work by Stilson and Smith [4], who
considered virtual analog modeling of the Moog VCF in the linearised case, and
Huovilainen, who examined the nonlinear case [5].

Consider the nonlinear Moog voltage-controlled filter system under non-forced,
autonomous conditions:

$\dot{\bf x}={\bf f}\left({\bf x}\right)\qquad\qquad{\rm
where}\qquad\qquad{\bf x}=\begin{bmatrix}x_{1}\\\ x_{2}\\\ x_{3}\\\ x_{4}\\\
\end{bmatrix}\qquad\qquad{\bf f}\left({\bf
x}\right)=\omega_{0}\begin{bmatrix}-\tanh(x_{1})-\tanh(\alpha^{4}x_{4})\\\
-\tanh(x_{2})+\tanh(x_{1})\\\ -\tanh(x_{3})+\tanh(x_{2})\\\
-\tanh(x_{4})+\tanh(x_{3})\\\ \end{bmatrix}$ (1)

for constant $\omega_{0}>0$, which serves as a cutoff frequency, and
$\alpha=\sqrt{2}r^{1/4}$, for constant $r$ with $0\leq r\leq 1$, which
controls the resonance of the filter. Here, ${\bf x}$ is the state,
representing the voltages (nondimensionalized) across a series of four
capacitors. ${\bf x}$ is a function of time $t$, and $\dot{\bf x}$ indicates
the ordinary time derivative of ${\bf x}$.

What we would like to find for this system is a Lyapunov function $V({\bf
x})$, satisfying

  * •

(C1) $V({\bf x})=0\qquad$ for $\qquad{\bf x}={\bf 0}$

  * •

(C2) $V({\bf x})>0\qquad$ for $\qquad{\bf x}\neq{\bf 0}$

  * •

(C3) $\dot{V}\leq 0$ for $\qquad{\bf x}\neq{\bf 0}$

The third condition above can be strengthened to asymptotic stability using a
strict inequality:

  * •

(C3)’ $\dot{V}<0$ for $\qquad{\bf x}\neq{\bf 0}$ (asymptotic stability)

If, furthermore, we have the condition of radial unboundedness, or

  * •

(C4) $V({\bf x})\rightarrow\infty\qquad$ as $\qquad|{\bf x}|\rightarrow\infty$
(global asymptotic stability)

we then have a condition for global asymptotic stability.

We also require that the Lyapunov function satisfy these properties for the
whole range of parameters $\omega_{0}>0$ and $0\leq r\leq 1$, though in this
short report, we do not assume time-varying behaviour.

##  2 Linearised System

Under linearisation, the Moog VCF system reduces to

$\dot{\bf x}={\bf Ax}\qquad\qquad{\rm where}\qquad\qquad{\bf
A}=\omega_{0}\begin{bmatrix}-1&0&0&-\alpha^{4}\\\ 1&-1&0&0\\\ 0&1&-1&0\\\
0&0&1&-1\\\ \end{bmatrix}$ (2)

For the linearised system, stability is entirely determined by the eigenvalues
$\\{\lambda_{{\bf A}}\\}$ of ${\bf A}$, which are

$\\{\lambda_{{\bf
A}}\\}=-\omega_{0}+\omega_{0}\alpha\\{e^{j\pi/4},\,e^{3j\pi/4},\,e^{5j\pi/4},\,e^{7j\pi/4}\\}$
(3)

Under the condition that $0\leq r\leq 1$, we then have
$0\leq\alpha\leq\sqrt{2}$, and it is true that

${\rm max}\left({\rm Re}\\{\lambda_{{\bf A}}\\}\right)\leq 0$ (4)

and thus all solutions are exponentially non-increasing for the full range of
parameters $0\leq r\leq 1$. We would like to determine a Lyapunov function
which shows the stability of this system under this same range of parameters.
(NB: In this case of a linear system, we could simply diagonalise the system
above, but we are interested in extensions to the nonlinear case.)

The natural choice of a Lyapunov function here is

$V({\bf x})=\frac{1}{2}{\bf x}^{T}{\bf x}$ (5)

$V$ corresponds to the physical stored energy across the four capacitors in
the Moog VCF, and obviously satisfies conditions C1, C2 and C4. It remains to
check conditions C3 and C3’. We have:

$\dot{V}=\dot{{\bf x}}^{T}{\bf x}={\bf
x}^{T}\underbrace{\left(\frac{1}{2}\left({\bf A}+{\bf
A}^{T}\right)\right)}_{{\bf A}_{s}}{\bf x}$ (6)

However, the matrix ${\bf A}_{s}$ is negative definite only for $0\leq
r<5/12$, for which it satisfies condition C3’ (and the weaker condition C3 if
$r=5/12$ is included). This choice of Lyapunov function does not recover the
full stable parameter range for the linearized Moog VCF.

A better choice can be achieved through scaling of the state variables.
Consider the new variable ${\bf w}$, defined as

${\bf w}={\bf D}{\bf x}\qquad{\rm where}\qquad{\bf
D}=\begin{bmatrix}1&0&0&0\\\ 0&\alpha&0&0\\\ 0&0&\alpha^{2}&0\\\
0&0&0&\alpha^{3}\\\ \end{bmatrix}$ (7)

(Here we leave aside the special case of $\alpha=0$, and assume
$0<\alpha\leq\sqrt{2}$.) Now we have

$\dot{{\bf w}}={\bf B}{\bf w}\qquad{\rm where}\qquad{\bf
B}=\omega_{0}\begin{bmatrix}-1&0&0&-\alpha\\\ \alpha&-1&0&0\\\
0&\alpha&-1&0\\\ 0&0&\alpha&-1\\\ \end{bmatrix}$ (8)

Now choose

$V({\bf w})=\frac{1}{2}{\bf w}^{T}{\bf w}=\frac{1}{2}{\bf x}^{T}{\bf
D}^{2}{\bf x}$ (9)

which again satisfies C1, C2 and C4. For conditions C3 and C3’,

$\dot{V}=\dot{{\bf w}}^{T}{\bf w}={\bf
w}^{T}\underbrace{\left(\frac{1}{2}\left({\bf B}+{\bf
B}^{T}\right)\right)}_{{\bf B}_{s}}{\bf w}$ (10)

Now, ${\bf B}_{s}$ is negative definite over the full range of parameters
$0<r<1$ and thus satisfies C3’ over this range (and the weaker condition C3 if
$r=1$ is included).

##  3 Nonlinear Moog VCF

Consider now the Moog VCF system as defined in (1). Define now a general
diagonal scaling of the form

${\bf w}={\bf D}{\bf x}\qquad{\rm where}\qquad{\bf
D}=\begin{bmatrix}1&0&0&0\\\ 0&d&0&0\\\ 0&0&d^{2}&0\\\ 0&0&0&d^{3}\\\
\end{bmatrix}$ (11)

for some real parameter $d>0$. (Earlier, in the linear case, we set
$d=\alpha$, but here we leave $d$ unspecified for the moment.) In terms of
${\bf w}$, the system may now be written as

$\dot{{\bf
w}}=\omega_{0}\begin{bmatrix}-\tanh(w_{1})-\tanh(\alpha^{4}w_{4}/d^{3})\\\
-d\tanh(w_{2}/d)+d\tanh(w_{1})\\\
-d^{2}\tanh(w_{3}/d^{2})+d^{2}\tanh(w_{2}/d)\\\
-d^{3}\tanh(w_{4}/d^{3})+d^{3}\tanh(w_{3}/d^{2})\\\
\end{bmatrix}=\omega_{0}\underbrace{\begin{bmatrix}-1&0&0&-d\\\ d&-1&0&0\\\
0&d&-1&0\\\ 0&0&d&-g\\\ \end{bmatrix}}_{{\bf
Q}}\underbrace{\begin{bmatrix}\tanh(w_{1})\\\ d\tanh(w_{2}/d)\\\
d^{2}\tanh(w_{3}/d^{2})\\\ \frac{1}{d}\tanh(\alpha^{4}w_{4}/d^{3})\\\
\end{bmatrix}}_{{\bf z}}$ (12)

where note ${\bf z}$ and ${\bf Q}$ as defined above, and where

$g(w_{4})=d^{4}\frac{\tanh(w_{4}/d^{3})}{\tanh(\alpha^{4}w_{4}/d^{3})}$ (13)

Notice the $g$ takes on values between $d^{4}$ and $d^{4}/\alpha^{4}$.

###  3.1 Lyapunov Function

Now, define a candidate Lyapunov function for the nonlinear Moog VCF system as

$V({\bf
w})=\ln(\cosh(w_{1}))+d^{2}\ln(\cosh(w_{2}/d))+d^{4}\ln(\cosh(w_{3}/d^{2}))+\frac{d^{2}}{\alpha^{4}}\ln(\cosh(\alpha^{4}w_{4}/d^{3}))$
(14)

(Recall that we have specified $d>0$ and also have restricted $\alpha$ such
that $0<\alpha\leq\sqrt{2}$.) $V$ again satisfies conditions C1, C2 and C4 for
any choice of $d>0$. For conditions C3 and C3’, we have

$\dot{V}=(\nabla V)^{T}\dot{{\bf w}}={\bf z}^{T}\dot{{\bf w}}$ (15)

using the gradient $\nabla V$ of $V$ with respect to ${\bf w}$, and where
${\bf z}$ is as defined in (12). But, using the definition of ${\bf Q}$ from
(12), we then have

$\dot{V}=\omega_{0}{\bf z}^{T}{\bf Q}{\bf z}=\omega_{0}{\bf
z}^{T}\underbrace{\left(\frac{1}{2}\left({\bf Q}+{\bf
Q}^{T}\right)\right)}_{{\bf Q}_{s}}{\bf z}$ (16)

It then remains to determine conditions under which ${\bf Q}_{s}$ is negative
definite. To this end, rewrite ${\bf Q}_{s}$ as

${\bf Q}_{s}=-{\bf I}_{4}+\frac{d}{2}\underbrace{\begin{bmatrix}0&1&0&-1\\\
1&0&1&0\\\ 0&1&0&1\\\ -1&0&1&f\\\ \end{bmatrix}}_{{\bf G}}$ (17)

in terms of a matrix ${\bf G}$, which is dependent on $f$, defined in terms of
$g$ as

$f=\frac{2}{d}\left(1-g\right)$ (18)

The maximal positive eigenvalue of ${\bf G}$ may be written as

$\lambda_{{\bf
G},max}=\max\left(\sqrt{2},\frac{1}{2}\left(f+\sqrt{f^{2}+8}\right)\right)$
(19)

To satisfy condition C3’, we want to find a functional form $d(\alpha)$ such
that

$\lambda_{{\bf G},max}<2/d\qquad{\rm for\,all}\qquad\alpha\qquad{\rm
with}\qquad 0<\alpha\leq\sqrt{2}$ (20)

It is useful at this point to consider two cases for $\alpha$:

#### Case I: $\alpha\leq 1$

In this case, choose $d=1$. Thus $g$ takes on values between $1/\alpha^{4}$
and $1$, so $g\geq 1$. This implies that $f\leq 0$, and thus $\lambda_{{\bf
G},max}=\sqrt{2}$, and thus (20) is satisfied.

#### Case II: $1<\alpha\leq\sqrt{2}$

In this case, we can choose $d=\alpha$. Thus $g$ takes on values between $1$
and $\alpha^{4}$, so $g\geq 1$. Again this implies that $f\leq 0$, and thus
$\lambda_{{\bf G},max}=\sqrt{2}$, and furthermore that

$\lambda_{{\bf G},max}=\sqrt{2}\leq 2/\alpha$ (21)

and again (20) is satisfied.

Thus, the choice

$d=\max(1,\alpha)$ (22)

leads to a Lyapunov function (14) that recovers the full range of values of
the parameter $\alpha$:

$\alpha\leq\sqrt{2}\qquad\rightarrow\qquad r\leq 1$ (23)

## References

  * [1] R. Moog.  A voltage-controlled low-pass high-pass filter for audio signal processing.  In 17th AES Convention, New York, Oct. 1965.
  * [2] T. Hélie.  On the use of Volterra series for real-time simulations of weakly nonlinear analog audio devices: Application to the Moog ladder filter.  In Proceedings of the 9th International Digital Audio Effects Conference, pages 7–12, Montreal, Canada, September 2006.
  * [3] T. Hélie.  Lyapunov stability analysis of the moog ladder filter and dissipativity aspects in numerical solutions.  In Proceedings of the 14th International Digital Audio Effects Conference, pages 45–52, Paris, France, September 2011.
  * [4] T. Stilson and J. O. Smith III.  Analyzing the Moog VCF with considerations for digitalimplementation.  In Proceedings of the International Computer Music Conference, pages 398–401, Hong Kong, August 1996.
  * [5] A. Huovilainen.  Nonlinear digital implementation of the moog ladder filter.  In Proceedings of the 7th International Conference on Digital Audio Effects, pages 61–64, 2004.