There are 9 plots, 5 from the paper pages itself, 2 simply graphic black versions of the same first two plots, and 2 more additional auxiliary plots (adaptive bin data cited in the paper with log-scale magnitude-like representation and sinusoid plus the visually-persuasive plot of LAT+optical superposed data).


Figure 1: LAT $\gamma$-ray lightcurves of \pg\ over $\sim$ 6.9 years, from 2008 August 4 to 2015 July 19. The lightcurve above 1 GeV is shown with a constant 45-day binning (top panel); two light curves above 100 MeV are shown, with 45- and 20-day binning (middle and bottom panels).


Figure 2: Multifrequency lightcurves of \pg\ at X-ray, optical and radio bands. Top panel: \textit{Swift}-XRT integrated flux (0.3-2.0 keV).Central panel:  optical flux density from Tuorla (R filter, black filled circle points), Catalina CSS (V filter rescaled, blue filled squared points), KAIT (V filter rescaled, red filled diamond points), \textit{Swift}-UVOT (V filter rescaled, green filled circle points). Dotted line:  LAT flux ($E > 100$ MeV) with time bins of 20 days, scaled and y-shifted for comparison.
Bottom panel: 15 GHz flux density from OVRO 40m (black filled circle points) and parsec-scale 15 GHz flux density by VLBA (MOJAVE program, yellow diamond filled points).


Figure 3: Left top panel: pulse shape (epoch-folded) $\gamma$-ray ($E>100$ MeV) flux lightcurve at the 2.18 year period (two cycles shown).
Left bottom panels: 2D plane contour plot of the CWT power spectrum (scalogram) of the $\gamma$-ray lightcurve, using a Morlet mother function (filled color contour).
This epoch-localized technique shows that most of the power is found at a timescale $ \sim 2$ years (clear/white horizontal strip) along all the MJD epochs, indicating a dominant periodic modulation.
The side panel to this is the 1D smoothed, all-epoch averaged, spectrum of the CWT scalogram showing a signal power peak in agreement with the 2.18-year value, also showing the LSP.
In this plot also the LSP is represented, with a peak in agreement with the average-CWT peak and the epoch-folding analysis. The statistical significance of the signal peak is evaluated against Brownian red (low frequencies) noise, through Monte Carlo first-order autoregressive process simulation.
Dashed lines depict increasing levels of confidence against red-noise calculated with Monte Carlo simulation.
The $\gamma$-ray signal peak is above the 99\% confidence contour level  ($<1$\% chance probability of being spurious).
Right top panel: pulse shape from epoch folding of the optical flux lightcurve at the 2.18 year period (two cycles shown).
Right bottom panels: the same CWT and LSP diagrams for the optical lightcurve.
The optical signal peak is above the 95\% confidence contour level.


Figure 4: Power Density Spectrum (PDS) of the LAT, $0.1-300$ GeV, count rate (photons cm$^{-2}$ s$^{-1}$)
lightcurve of \pg\ from a 3$^{\circ}$ exposure-weighted aperture photometry technique with 600-second time bins.


Figure 5: Discrete cross-correlation plots from the approach with PDS model measured from the lightcurve data \citep{maxmoerbeck14}. In each plot the black dots are the DCCF estimates, and the red, orange and green lines are the  1$\sigma$, 2$\sigma$ and 3$\sigma$ significance levels respectively. Top panel: DCCF between the radio 15GHz and $\gamma$-ray (20-day time bins) lightcurve.
Central panel: DCCF between the unbinned optical lightcurve and $\gamma$-ray (20-day time bins) lightcurve.
Bottom panel: DCCF between the 20-day rebinned optical lightcurve and $\gamma$-ray (20-day time bins) lightcurve.
The oscillating shape of the significance contours for this case is due to the number of samples in each bin.


Figure 1 black version.


Figure 2 black version.
 
 
Figure 6 auxiliary: composite optical-gamma-ray flux light curve of PG 1553+113 highlighting the 2.18-year modulation. Optical data are all the cumulated Tuorla, Catalina CSS, KAIT and UVOT data points (gray open symbols) reported in Figure 2. Fermi LAT data (E > 100 MeV, 20-day bins, re-scaled flux) are represented by red filled symbols.
 
 
Figure 7 auxiliary: gamma-ray flux light curve (E > 100 MeV) of PG 1553+113 in 20-day bin (blue open circles, the same data of Figure 1, likelihood technique) superposed to the LAT adaptive-binning flux light curve (black line). The two light curves are here presented in a logarithmic scale to mimic the common plots for the optical magnitude light curves of variable stars. A strictly-periodic sinusoidal waveform of period 2.18 years is superposed to the data.