text:
\section{Introduction}
The discovery of the energetic, stripped-envelope core-collapse SN~1998bw, which was temporally and spatially coincident with GRB~980425 (Galama et al. 1998; Patat et al. 2001), provided the first empirical evidence that GRBs arise from massive stars.  Since 1998 there have been several more spectroscopic (e.g. Pian et al. 2006; Sparre et al. 2011; Bufano et al. 2012; Xu et al. 2013; Schulze et al. 2014) and photometric observations (e.g. Stanek et al. 2005; Bersier et al. 2006; Cano et al. 2011B) of supernovae (SNe) accompanying GRBs (GRB-SNe).  GRB-SNe are very energetic, possessing, on average, ten times more kinetic energy than SNe Ibc not associated with GRBs (Iwamoto et al. 1998).  Every single GRB-SNe that has been spectroscopically classified is of type Ic, implying that their outer layers of hydrogen and helium have been lost prior to explosion.  Moreover, the blueshifted line velocities in their spectra are much larger than those measured for non-GRB SNe Ibc, approaching velocities more than 10\% the speed of light, and which display broad-lined (BL) absorption lines.  Such rapid line velocities and broad-lined profiles have been seen in the spectra of SNe Ic-BL not associated with GRBs, which are also more energetic and have larger ejecta velocities than SNe Ibc, though not as extreme as GRB-SNe (Cano 2013).

The exact nature of the central engine in GRB-SNe is uncertain, though theory suggests it may being either an accreting black hole, or a rapidly rotating neutron star that possesses a very strong magnetic field.  In these models (e.g. Woosley 1993; MacFadyen et al. 1999; MacFadyen et al. 2001), a stellar black hole is created during the core-collapse of a rapidly rotating massive star (e.g. Yoon \& Langer 2005; Woosley \& Heger 2006).  The angular momentum present in the core at the time of collapse leads to the formation of a torus of material that accretes onto the BH, ultimately producing an ultra-relativistic, bi-polar jet along the rotation axis.  Once the jets pierce through the star, they interact with themselves producing the initial burst of $\gamma$-rays, and later they collide with the surrounding medium, producing the longer-lasting X-ray/optical/radio afterglow (AG).  It has been suggested by Bromberg et al. (2011) that low-luminosity GRBs such as GRB~060218 and GRB~100316D arise from jets that just barely breakout from the stellar surface, and therefore emit less $\gamma$-ray radiation.

Adding diversity to this scenario is the existence of type Ic-BL SN~2009bb (Soderberg et al. 2010; Pignata et al. 2011).  To date, SN~2009bb is one of only two SN Ic-BL (the other being SN~2012ap; Margutti et al. 2014; Milisavljevic et al. 2014) not associated with a GRB for which there is compelling observational evidence that they both possessed a central engine that contributed to the explosion of their stellar progenitors.  Luminous radio detections of SN~2009bb were attributed by Soderberg et al. (2010) as arising from a substantially relativistic outflow powered by a central engine.  The kinetic energy contained within the relativistic outflow of SN~2009bb was demonstrated to be similar to that of low-luminosity GRBs 060218 and 100316D (Margutti et al. 2013; Margutti et al. 2014), which indirectly suggests that SN~2009bb may have a similar degree of collimation/asymmetry as low-luminosity GRB-SNe.  In the above theoretical scenario, SN~2009bb arises from a jet that failed to breakout of the star, possibly because of a less energetic central engine.

Another open question is whether GRBs and/or GRB-SNe can be used as standardizable candles.    While a few GRB-SNe have occurred in the local universe, with the relative brightness of the event allowing a direct analysis of their observational and physical properties to be performed, isolating the SN contribution from the optical and near-infrared (NIR) observations of fainter events is more complicated.  Precisely ascertaining the observational properties of the accompanying SNe requires a careful decomposition of the photometric observations.  The degree of accuracy obtained in decomposition methods depends strongly on precisely knowing the total amount of extinction along the line of sight (both foreground and rest-frame), where determining the latter requires assumptions that complicate the analysis.  Furthermore, the contribution from the host galaxy needs to be quantified and removed, either by image-subtraction techniques, or by mathematically subtracting the host flux from the observations.  Finally, in order to isolate and model the flux from the SN, the temporal behaviour of the optical AG, which arises from the initial GRB event, needs to be precisely determined.  If the AG is improperly modelled, the AG contribution can be under- or over-estimated, which affects the resultant SN properties.  Therefore, in order for the SN to be modelled accurately, each one of these hurdles needs to be carefully overcome.

Several studies (Zeh et al. 2004; Ferrero et al. 2006; Cano et al. 2011A; Th{\"o}ne et al. 2011; Cano 2013) have attempted to determine the luminosity ($k$) and stretch ($s$) factor of GRB-SNe relative to a template SN, which historically has been SN~1998bw.  These parameters are analogous to the absolute peak SN magnitude and $\Delta m_{15}$ (the amount the light curve fades from peak light to fifteen days later) used in the luminosity--decline relation for SNe Ia (Phillips 1993).  Initial studies such as those by Ferrero et al. (2006) and Cano et al. (2011A) noted that no obvious correlation exist in the observer-frame $R$-band $k$- and $s$-factors of several GRB-SNe.  However, as the GRB-SNe in their respective samples occurred over a broad range of redshifts (0.0085$\le z \le1.006$), the $R$-band light curves (LCs) in their respective analyses were all probing different rest-frame wavelengths, and thus different regions of their spectral energy distributions (SEDs).  Encouragingly, a recent analysis by Schulze et al. (2014) looked at the \textit{approximate} rest-frame $V$-band stretch and luminosity factors of several GRB-SNe, finding that a rough correlation was present.  An additional analysis by Lyman et al. (2014) also hinted at a correlation between the peak absolute magnitude and $\Delta m_{15}$ of several nearby GRB-SNe and SN~2009bb.  Excitingly, a very recent paper by Li \& Hjorth (2014) have presented a luminosity--decline correlation between the K-corrected\footnote{Not to be confused with the $k$ parameter used here to denote the brightness of a given GRB-SNe in a given filter, but it is rather the correction that arises from observing events occurring at vast cosmological distances, where the observer-frame light arises from bluer rest-frame light that has been redshifted as it traverses through an expanding universe.}, $V$-band, peak absolute magnitude and decay rates from peak of eight GRB-SNe.  Reassuringly, their results fully support those presented in this paper.

In this present work I performed a detailed analysis to precisely account for the varying redshifts of GRB-SNe in order to determine the \textit{exact} rest-frame ($UBVRI$) stretch and luminosity factors of several GRB-SNe.  The latest values of the observer-frame (foreground) extinction due to the Milky Way (MW) from Schlafly \& Finkbeiner (2011) have been used, as well as values for the rest-frame extinction derived by several works (see the Appendix for a detailed discussion of each GRB-SN).  Importantly, this analysis is an improvement over previous works (e.g. Cano et al. 2011A; Cano 2013) by calculating the \textit{redshift-corrected} rest-frame extinction using the value of $A_{V,\rm rest}$, and the MW/SMC/LMC templates of Pei (1992).  This is more precise than considering a total line-of-sight extinction for a given filter, for at larger redshifts, the extinction measured in a certain observer-frame filter will arise from bluer light in the rest-frame that will be correspondingly more extinguished.  A $\Lambda$CDM cosmology constrained by Planck (Planck Collaboration et al. 2013) has been employed here, where $H_{0} = 67.3$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_{M} = 0.315$, $\Omega_{\Lambda} = 0.685$.

\section{Procedure}

\subsection{Decomposing the optical light curves}

The method used involves decomposing the optical LCs to isolate flux coming from just the accompanying SN.  Previous papers (e.g. Zeh et al. 2004; Ferrero et al. 2006; Cano et al. 2011A; Th{\"o}ne et al. 2011; Cano 2013) have used similar methods to determine the brightness ($k$) and width/stretch factor ($s$) in a particular filter relative to a template SN, which historically has been SN~1998bw due to its well-sampled LC and well-determined temporal evolution.  A step-by-step visual diagram of the decomposition method is displayed in Fig. \ref{fig:090618}.

For each GRB-SN in the sample, I collected all of the available data published in the literature (Table \ref{table:refs}).  While the amount of SN-bumps in optical and NIR LCs is more than 20 (Cano 2013), the analysis has been limited to those events where: (1) observations were in at least two filters in order to interpolate between filters to get to the precise to rest-frame wavelength), (2) observations of the host had been obtained, (3) knowledge of the rest-frame extinction had been determined, and (4) the LC is well sampled enough that the AG and SN temporal behaviour can be accurately and, importantly, unambiguously modelled.  Using this criteria, eight GRB-SNe were analysed, not including the template SN~1998bw.  Once the data were gathered, the magnitudes were corrected for foreground extinction and then converted into monochromatic fluxes, using the effective wavelengths and flux zeropoints in Fukugita et al. (1995).

A vital part of the decomposition is quantifying the total extinction along the line of sight.  Observer-frame (foreground) extinction due to various sight-lines through the Milky Way was corrected for using the dust maps of Schlafly \& Finkbeiner (2011).  The rest-frame extinction is more complex to determine.  One way is to obtain a spectrum of the host galaxy after the GRB and SN have faded away, and use either absorption or emission line diagnostics to estimate the host-galaxy extinction.  The draw-back of this method is that for GRBs at higher redshift that cannot be resolved in their host galaxy, one has to make the assumption that the integrated dust extinction of the entire galaxy is a suitable proxy for the dust-extinction local to the GRB event itself, which for some GRB events has proved to be a poor proxy (e.g. Kr{\"u}hler et al. 2011).  Another method is to model early observations of the AG at optical, IR and X-ray frequencies, and fit the data with extinction curves based on the SMC, LMC and MW (e.g. Kann et al. 2006; Guidorzi et al. 2009; Kann et al. 2010; Cano et al. 2011B; Olivares et al. 2012).  The rest-frame extinction for the majority of the GRB-SNe in the sample has been determined using the latter method.  Note that the rest-frame extinction is only corrected for after (1) the observer-frame host contribution has been removed, and (2) the interpolation described below was completed.

For each GRB-SN, the light is attributed as coming from three different sources: (1) the AG, which is associated with the GRB event, (2) the SN, and (3) the constant source of flux coming from the host galaxy.  The host contribution was determined and removed either by image subtraction (e.g. GRB 011121; Bloom et al. 2002), or by observing the host galaxy at late times after the AG and SN have faded away, and then mathematically subtracting the host flux from the earlier epochs (e.g. GRB 130831A; Cano et al. 2014).  Note that the host galaxy contribution needs to be removed at this stage, because once the SED interpolation is performed, we do not have knowledge of the brightness of the host galaxy at the redshift-corrected, observer-frame wavelength to later remove.  After removing the host contribution, the remaining data is attributed to coming from just the AG and SN.

At this point, observer-frame SEDs at each epoch of contemporaneous optical data for each GRB-SN are created.  For example, for GRB 090618 (redshift $z=0.54$), I have used observer-frame $R_{c}$ and $i$ data from Cano et al. (2011B), which using the effective wavelengths of each filter from Fukugita et al. (1995) corresponds to observer-frame wavelengths of 6588 and 7706 \AA $ $ respectively.  Rest-frame $B$-band is 4448 \AA, which corresponds to observer-frame 4448($1+z$) $=$ 6850 \AA.  Therefore at each epoch of host-subtracted data,  the SED has been interpolated to this wavelength to extract the flux at redshifted rest-frame $B$-band.

Once the redshift-corrected LC is created, which has now been corrected for rest-frame extinction, it needs to be modelled to determine the temporal behaviour of the AG, as well as the properties of the accompanying SN.  Standard GRB theory suggests that the light powering the AG is synchrotron in origin, and therefore follows a power-law behaviour in both time and frequency ($f_{\nu} \propto (t - t_{0})^{-\alpha}\nu^{-\beta}$, where the respective decay and energy spectral indices are $\alpha$ and $\beta$, while $t_{0}$ is the trigger time of the GRB).  It is therefore very beneficial to have a wide time coverage of the entire event in order to accurately determine the rate at which the AG decays, especially at late times when the SN is the dominant source of flux.

Finally, in order to measure the luminosity and stretch factors of each accompanying SN, I used a redshifted (K-corrected) template of SN 1998bw, originally presented in Cano (2013).  The templates are created using a \textsc{c}-program that calculates the rest-frame filter LCs for each SN and interpolates in frequency and time over the original LCs of SN~1998bw (Galama et al. 1998; Patat et al. 2001), and then calculates the luminosity distance using a $\Lambda$CDM cosmology to rescale the flux density.  The program generates a synthetic LC at the desired observer-frame wavelength.

The synthetic LC is then fitted with a log-linear spline using \textsc{pyxplot}\footnote{http://pyxplot.org.uk}.  If we denote the log-linear spline as function $\Lambda(x)$, then to determine the stretch and luminosity factors of each SN relative to the template we use a new function of the form $\Sigma = k \times \Lambda(x/s)$.  This function is added to the AG function, and the data is fitted to simultaneously determine the decay rate of the AG at all times, as well as $k$ and $s$, all of which are free parameters.  It is seen that the LC of GRB~090618 has a ``break'' in it, where a steepening of the AG occurs at some break time ($T_{\rm B}$).  Thus, two additional free parameters are included ($T_{\rm B}$ and the decay constant after the break).  A summary of the best-fit parameters are presented in Table \ref{table:params}.

Throughout this analysis there are several sources of errors: host-galaxy subtraction, afterglow subtraction, uncertainties in the observer-frame and rest-frame extinction and the SED construction and interpolation.  Each one of these procedures/quantities introduces an error into the final calculation of $k$ and $s$.  The errors of the first four are photometric, and the errors have been added in quadrature.  A conservative estimation of the total error in $k$ and $s$, which includes an estimate of the error introduced by the SED interpolatiom, is of order 20\%, and is the fractional value adopted in Table \ref{table:params}.

A visual example of the error propagation can be seen in Fig. \ref{fig:090618}, where the final errorbars of the rest-frame $B$-band LC are larger than those of the observer-frame $R$- and $i$-band LCs as they also include the uncertainty associated with the rest-frame extinction.

\section{Results}

For every GRB-SNe in this sample, I interpolated the optical/NIR spectral energy distributions (SEDs) to exact rest-frame frequencies, and then performed the LC decomposition method to determine $k$ and $s$ at every possible rest-frame filter.  The resultant value of $k$ and $s$ in each rest-frame filter provides an independent measure of the shape of the SN relative to the template SN.  Note that this approach can be considered to be model-dependent as it relies on using a template SN, rather than focusing on the SN properties in and of themselves.  As such, I have assumed that each GRB-SN evolves in a similar fashion to the template SN, corrected only for the unique intrinsic brightness and temporal evolution of each GRB-SN, but with similar color evolution.  It has been shown by e.g. Lyman et al. (2014) that GRB-SNe, and all SNe for that matter, will possess a unique colour evolution.  However, it should be noted that a GRB-SN with $k$ and $s$ values that are not identical do account for their unique colour evolution relative to the template SN.  For example, a fiducial event with $s_{R}=0.8$ and $s_{I}=0.9$ implies that the SN gets redder over time more rapidly than an event with $s_{R}=s_{I}=0.9$ (in the former event, $R-I$ will increase more rapidly after peak light as the SN fades more rapidly in $R$ relative to $I$).  As such, the method used here does provide a way to account for the unique colour evolution of each GRB-SN.

With this caveat in mind, the analysis of Li \& Hjorth (2014) can be considered to be less model dependent in this respect, though it should be considered that their analysis still relies on using the SED of SN~1998bw as a template when calculating the K-correction of those events in their sample where their constructed SED is incomplete.

Next, using \textsc{python} programs that employ \textsc{numpy} and \textsc{scipy} algorithms, I performed a non-linear least squares analysis on the combined dataset to find the best-fitting values of the slope ($m$) and y-intercept ($b$) to a straight line that was fitted to the data.  The fit was performed with the \textsc{optimization.curve$\_$fit} routine from \textsc{scipy}, which uses the Levenberg-Marquardt algorithm to fit the line to the data. The analysis was designed to (necessarily) consider the error in both $k$ and $s$.  To obtain an estimate of the error in $m$ and $b$, I performed a Monte Carlo analysis, where a new dataset is created from the original by choosing a new value from a Gaussian that is centered on each original datapoint, and whose standard deviation in each direction is equal to the error in that variable.  The simulation was performed 10,000 times, each time fitting the straight line to the new dataset.  The estimated errors of $m$ and $b$ arise from their resultant distributions, where the quoted 1-$\sigma$ errors are the standard deviation of each distribution.  The best-fitting values are $m=1.60\pm0.20$ and $b=-0.31\pm0.15$.  Interestingly, this result implies that GRB-SNe do not follow a 1:1 ratio of $k$:$s$, rather $k$ increases more rapidly with larger $s$.  The results are displayed in Figure \ref{fig:GRB-SNe}.

It is confirmed that a statistically significant correlation is present in the data.  The Pearson's correlation coefficient of the sample ($N=9$) is $r=0.936$, while the two-point probability of a chance correlation is $p=4.2\times10^{-9}$.  It is worth noting that every single nearby, spectroscopically identified GRB-SN obeys this relation, including the ``Big Five'' GRB-SNe (GRBs 980425, 030329, 031203, 060218 and 100316D) and GRB~120422A.  The notable absentee is recent GRB 130427A, which is associated with SN~2013cq (Xu et al. 2013).  The exclusion of this event is dictated by its complex AG behaviour and uncertain host galaxy contribution (Perley et al. 2014), which, like GRB 030329, severely hinders our ability to precisely decompose the optical and NIR observations, and will most likely require more advanced decomposition methods such as those used by Deng et al. (2005) on the spectra obtained of GRB 030329.

Motivated by the observation that SN~2009bb is also a relativistic, type Ic-BL SN, I included its $BVRI$ stretch and luminosity factors from Cano (2013) with the GRB-SNe sample.  I have adopted an uncertainty of 20\% in the values of $k$ and $s$ for SN~2009bb in this analysis, which are larger than those originally presented in Cano (2013).  I then refit the data using the non-linear least squares and Monte Carlo analysis (Figure \ref{fig:SNe_GRB_IcBL}).  The Pearson's correlation coefficient for this new sample is $r=0.915$, with a chance probability of $p=9.7\times10^{-10}$.  Moreover, the values of $m=1.53\pm0.19$ and $b=-0.31\pm0.14$ are in good agreement with those determined from just the GRB-SN sample.  This result suggests that perhaps all engine-driven SNe are potential standardizable candles, though this conclusion is understandably tentative as it is based only on a single non-GRB-SN event.

Next I considered the sample of non-GRB SNe Ibc originally presented in Cano (2013), where I determined $k$ and $s$ relative to SN~1998bw, to search for a similar correlation.  Only those events where the total extinction along the line of sight was known, and for which the host galaxy contribution was also removed, were included in the analysis.  Note that the stretch factors of the SNe Ibc are not corrected for redshift effects, however the bulk of the SNe Ibc considered here have redshifts of order $z\sim0.01$.  The resultant effect on the transformation between observer-frame to rest-frame for these very nearby events is of order a few percent, and thus does not fundamentally affect the results.  Plotted in Figure \ref{fig:SNe_Ibc} are the $k$- and $s$-factors of different SNe Ibc sub-types (Ib, Ic and Ic-BL), along with the best-fitting line, and 2-$\sigma$ uncertainty region found for the GRB-SNe sample.  Also displayed is the Pearson's correlation coefficient for each sample, as well as the combined SNe Ibc sample.  The low value of $r$ for each SN subtype highlights their unsuitability as standardizable candles.

\section{Discussion}

It has been noted that likely all core-collapse SNe possess a degree of asphericity (Wang \& Wheeler 2008), where the asphericity arises from, and is driven by, Rayleigh-Taylor instabilities surrounding the collapsing core.  This is in contrast to GRB-SNe, where asphericity also arises from the relativistic jet that is launched by the central engine.  Observationally, asphericity in SNe Ibc has been determined via different methods, such as modelling of optically thin nebular emission lines (e.g. Mazzali et al. 2005; Maeda et al. 2008; Taubenberger et al. 2009), where the emission lines of elements such as oxygen trace the overall velocity distribution of the ejecta. Another method is comparing the results of hydrodynamical explosion \& radiative transfer models with observations, in order to determine the ejecta geometry.  For example, it was shown that a large excess of luminosity along the polar/rotation axis relative to the equatorial axis (ratio of 4:1) was needed to explain the early bolometric LC of SN~1998bw (Maeda et al. 2006).  It is very probable that all GRB-SNe are highly aspherical (Woosley \& Bloom 2006; Maeda et al. 2006).  However, because we detect the initial $\gamma$-ray emission, as observers we must be near to the jet-axis, or the high-energy emission would not be observed.  The subsequent SN therefore must also be viewed close to the jet axis, meaning that we are observing each event at approximately the same viewing angle.  Conversely, because the non-GRB SNe Ibc are viewed at all possible viewing angles, and with each possessing a different degree of asphericity, it is not unexpected that no relation exists between the luminosity and stretch factors of these SNe Ibc relative to the template SN.  In the case of SN~2009bb, which was observed to be an engine-driven SN with a small amount of material moving at relativistic velocities similar to those of the low-energy GRB-SNe, the results here reinforce the supposition made by Pignata et al. (2011), that SN~2009bb was likely observed close to the polar/elongated axis.  It is tempting then to conclude that the possible use of GRB-SNe as standardizable candles arises, at least in part, from them possessing a similar ejecta geometry that is observed at roughly the same viewing angle.

Moreover, as the LCs of type I SNe (including GRB-SNe) are powered solely by the radioactive decay of $^{56}_{28}$Ni~$\rightarrow$~$^{56}_{27}$Co~$\rightarrow$ $^{56}_{26}$Fe (Arnett 1982; Woosley \& Weaver 1986), the fact that GRB-SNe have a distribution of brightnesses (i.e. $k$ values) results from each SN nucleosynthesizing a different amount of $^{56}$Ni.  If GRB-SNe are accurately described by the collapsar model, then the SN is expected to be driven by energy deposited by the jet as it disrupts the star, as well as by an energetic wind driven by viscous interactions in the accretion disk surrounding the newly formed compact object (MacFadyen et al. 1999).  Explosive nucleosynthesis of $^{56}$Ni occurs deep in the star as material displaced by the jet moves away at supersonic speeds, while further nucleosynthesis is driven by the wind from the accretion disk, with the latter possibly being the dominant producer of $^{56}$Ni for events such as SN~1998bw (MacFadyen et al. 2001).  A n\"aive conclusion to be drawn is that the correlation presented here suggests that a relation exists between the strength and energetics produced by the central engine and the resultant nucleosynthetic yields of $^{56}$Ni.  Moreover, the lack of correlation of $k$ and $s$ for the non-GRB SNe implies that the explosion and nucleosynthesis mechanism(s) are not correlated.

This n\"aive scenario does not account for the physical conditions of the SNe however.  In the context of SNe Ia, which are, of course, also standardizable candles, their LCs are also powered solely by the radioactive decay of nickel and cobalt, the amount of which determines the LC's peak brightness and width.  The width also depends on the photon diffusion time, which in turn depends on the physical distribution of the nickel in the ejecta, as well as the mean opacity of the ejecta.  In general, the opacity increases with increasing temperature and ionization (Woosley et al. 2007), thus implying that more nickel present in the ejecta leads to larger diffusion times.  This directly implies that fainter SNe Ia fade faster than brighter SNe Ia, thus satisfying the luminosity--decline relation (Phillips 1993).  This is not the only effect however, as the distribution of nickel in the ejecta also affects how the LC evolves, where nickel located further out has a faster bolometric LC decline.  Additionally, following maximum $B$-band light, SNe Ia colours are increasingly affected by the development of Fe II and Co II lines that blanket/suppress the blue $B$-band light.  Dimmer SNe are thus cooler, and the onset of Fe III $\rightarrow$ Fe II recombination occurs quicker than in brighter SNe Ia, resulting in a more rapid evolution to redder colours (Kasen \& Woosley 2007).  Therefore the faster $B$-band decline rate of dimmer SNe Ia reflects their faster ionization evolution, and provides additional clues as to why fainter SNe Ia fade more rapidly.  Thus, as the LCs of GRB-SNe are also powered by radioactive decay, the physics that govern SNe Ia also govern those of GRB-SNe, and may go some way to explaining why GRB-SNe are also standardizable candles.

\section{Conclusions}

In this paper I have shown that a statistically significant correlation exists between the brightness ($k$) and width/shape ($s$) of the LCs of GRB-SNe (the correlation coefficient being $r=0.932$) relative to a template SN (SN~1998bw).  This result highlights their potential use as standardizable candles.  Every single nearby, spectroscopically identified GRB-SN obeys this relation, including GRBs 980425, 030329, 031203, 060218, 100316D, 120422A and 130831A.  GRB~130427A has been omitted from this analysis due to its complicated LC structure, that will likely require more sophisticated decomposition such as that used by Deng et al. (2005) for GRB~030329.

When the relativistic, engine-driven type Ic-BL SN~2009bb (which is not associated with a GRB) was included in the analysis, it too was seen to follow this correlation ($r=0.916$ for the entire sample).  This highlights the possibility that perhaps all engine-driven SNe may eventually be used as standardizable candles, though this conjecture is very tentative as it is based on only a single event.

In this analysis, there are several sources of error that need to be considered.  These errors arise from the AG and host-galaxy subtractions, the uncertainties in the observer-frame and rest-frame extinction, and the SED interpolation.  In an attempt to properly account for the cumulative effect these will have on the resultant correlation between $k$ and $s$, I have adopted a conservative error of 20\% in the values of $k$ and $s$.  The adoption of this large error estimate is justified when compared with the original analysis of SNe Ia in Phillips (1993).  All nine SNe Ia in his sample were located at distances closer than SN~1998bw, which itself is the closest GRB-SN to date.  As such, the SNe Ia were all brighter than the GRB-SNe in the sample investigated here, and all had well-determined line-of-sight extinction values.  Moreover, no uncertainties were introduced from AG subtraction or SED construction and interpolation when determining the peak SN magnitudes and $\Delta m_{15}$ values.  This all implies that the luminosity--decline relation in Phillips (1993) is naturally better determined than the complicated decomposition technique employed here.

With one eye on the future, the modest size of the sample in this work highlights the need for additional observations of GRB-SNe, especially to obtain high-cadence and very importantly, multi-wavelength observations of each event.  The addition of further high-quality investigations of GRB-SNe will hopefully lead to a refinement and improvement of the scatter in the $k$--$s$ relationship, which in turn will allow for the possibility for their eventual use as a standardizable candle.  Once the redshift of a GRB has been determined, the resultant campaign to observe the accompanying SN should consider the rest-frame wavelengths and obtain observations to bracket as many rest-frame filters as possible.  Finally, but very importantly, in order to accurately determine the brightness and shape of the accompanying SN, it is necessary to have a good time coverage of the event at all epochs, as an incorrect AG model will give incorrect SN properties.

\acknowledgments

I am exceedingly grateful to Palli Jakobsson, Cristiano Guidorzi, David Alexander Kann and the anonymous referee for their comments and council on every aspect of this manuscript.  I happily acknowledge support by a Project Grant from the Icelandic Research Fund.


text:
\section{Analytical approach}
\label{sec:analytical_approach}

In this section we provide a detailed description of how to apply the framework, introduced in the main text, to approach gate design problems analytically.

In~\cref{subsec:perfect_state_transfer} we show how our framework can be applied to tackle generic perfect state transfer problems.
In~\cref{subsec:cnot_proof} we prove that a CNOT gate \textit{cannot} be implemented using only \emph{single-qubit} interactions.
While this is an obvious result, the calculation can be interesting to show, in a simple case, how the eigenvalue conditions given in~\cref{prop3} can be used to \emph{rule out} that a gate can be implemented with a specific set of interactions.
In~\cref{subsec:toffoli_main} we show how to use our framework to derive Hamiltonians implementing the Toffoli gate, using only one- and two-qubit interactions.
In~\cref{subsec:toffoli_posteriori_derivation} we show how one of the solutions obtained in~\cref{subsec:toffoli_main} could have been obtained via direct analytical reasoning, without any knowledge of the conditions given in the main text.
This provides some insight into how the solutions actually generate the Toffoli gate,
and illustrates the kind of ad-hoc non-trivial reasoning that, without the use of~\cref{prop123}, would have been necessary to find such solutions.

For notational convenience, we will in this section denote Pauli matrices with $X_i, Y_i, Z_i$, instead of $\sx{i}, \sy{i}$, and $\sz{i}$ as in the main text.

\subsection{Perfect state transfer}
\label{subsec:perfect_state_transfer}
A one-dimensional quantum walk is described by the Hamiltonian
$\HH_W =\sum_{k=1}^{N-1} J_k \ket{k}\bra{k{-}1} + B_k\ket k\bra k$,
where $N$ is the length of the lattice upon which the walk takes place, $J_k$ the transition rates between adjacent sites,
$B_k$ the local energies and
$\ket{k}$ defines the state where the ``walker'' is at the $k$-th site.
A quantum walk Hamiltonian $\HH_W$
admits perfect state transfer (PST) at time $t$, i.e. the initial state of the walker initially at site $k$ is perfectly retrieved at site $N-j+1$ if $e^{-it \HH_W} = \Xi$,
where $\Xi_{kj}=\delta_{k,N-j+1}$ is the reflection matrix.
Necessary and sufficient conditions for PST are well understood
\cite{yung2005perfect,kay2009a}: firstly, $\HH_W$ has to be ``mirror-symmetric'', that is
$[\HH_W,\Xi]=0$, and secondly the eigenvalues $\{E_k\}$ of $\HH_W$ should to satisfy the condition $e^{iE_kt}=(\pm 1)^k$.

We show now that finding the parameters $\{J_k, B_k\}$ for PST is a particular case of the Hamiltonian design problem for gate simulation.
Specifically, the conditions for PST can be
obtained from the construction in~\cref{prop123}. In a 1D quantum walk,
the ``physical'' couplings are the nearest neighbour interactions, but
$\HH_\Xi=i\log(\Xi)= \pi(\Xi-\mathds 1)/2$  is long range. Using $\HH_W$ as
$\HPhysLambda$, where $\lambda=\{J_k,B_k\}$, and defining $\HH' = \HH_W-\HH_\Xi$ following condition
\eqref{prop3}, where the physical interactions in $\HH_\Xi$ have been reabsorbed into the
definition of $\HH_W$, one finds that: (i)~condition \eqref{prop2} is equivalent to the mirror-symmetry request
$[\HH_W,\Xi]=0$; (ii) from conditions~\eqref{prop3} and from the definitions of $\mathcal{H}_W$ and $\mathcal{H}_\Xi$, one finds that the spectrum of
$\HH_W$ satisfies  $e^{iE_kt}=(\pm 1)^k$ --- indeed the eigenvalues of $\Xi$
are $\{0,\pi\}$, so $E_k=\pi(2n_k+1)$ for integer $n_k$.
It is straigthforward to extend the above proof for more general transfers, such as with
long-range interactions \cite{kay2006perfect}, or for perfect fractional revivals
\cite{banchi2015perfect,genest2016quantum}.

\subsection{Proof that CNOT needs 2-qubit terms}
\label{subsec:cnot_proof}
We here show how to use our framework to prove that a CNOT gate cannot be implemented using only one-qubit interactions.
While this result is trivial, it is nonetheless interesting to show how the framework can be used to obtain this kind of impossibility results.

The spectral decomposition of the CNOT reads
\begin{equation}
	\operatorname{CNOT} = Z_1^+ + Z_1^- X_2^+ - Z_1^- X_2^-,
	\label{eq:CNOT_definition}
\end{equation}
where we made a canonical choice for the basis of the three-fold degenerate eigenspace corresponding to the eigenvalue $+1$,
and defined $Z_i^\pm\equiv(1\pm Z_i)/2$, and similarly for $X_i^\pm$.
More explicitly, we are considering the following basis set of trace-1 projectors for the degenerate space: $\{Z_1^+ Z_2^+, Z_1^+ Z_2^-, Z_1^- X_2^+ \}$, whereas the fourth projector is bound to be $Z_1^- X_2^-$.
The corresponding principal Hamiltonian $\mathcal H_{\on{CNOT}}$, obtained by directly taking the logarithm of~\cref{eq:CNOT_definition}, is:
\begin{equation}
	\HH_{\on{CNOT}} =
	\pi Z_1^- X_2^-.
\end{equation}
Let us now consider what happens when the multivaluedness of the logarithm is taken into account, but no rotation of the degenerate eigenspace is performed.
Considering only the factors with two-qubit interactions, the following expression is found:
\begin{equation}
\begin{aligned}
	\HH_{\on{CNOT}} / 2\pi \sim \nu_{12} Z_1 Z_2 + 2\pi ( 1/2 + \nu_{43} ) Z_1 X_2,
\end{aligned}
\end{equation}
where here $\nu_{ij}\equiv \nu_i - \nu_j$, and $\nu_i\in\mathbb Z$ is the integer produced by application of the logarithm to the $i$-th projector.
Note how the $Z_1 X_2$ factor cannot be removed by \emph{any} choice of $\nu_i$,
which could be interpreted as a proof that two-qubit interaction terms are indeed necessary to implement the $\text{CNOT}$ gate.
This, however, does not in principle preclude the possibility that an appropriate rotation of the degenerate space allows to obtain a generator with only local terms.
To verify that this is not the case, we would have to consider a generic rotation $R$ of the degenerate space,
that is, an operator of the form
$R = \sum_{i,j=1}^3 r_{ij} \ketbra{+1_i}{+1_j}$,
with $\ket{+1_i}$ the $i$-th eigenvector in a fixed base of the degenerate space.
The problem is then that of finding a unitary $R$ and integers $\nu_i$ such that
\begin{equation*}
\begin{gathered}
	\nu_1 R (Z_1^+ Z_2^+) R^\dagger +
	\nu_2 R (Z_1^+ Z_2^-) R^\dagger +
	\nu_3 R (Z_1^- X_2^+) R^\dagger \\
	+ (\nu_4 + 1/2) Z_1^- X_2^-
\end{gathered}
\end{equation*}
does not contain 2-qubit interactions.
The solution of this problem is non-trivial,
mostly due to the many (9 in this case) parameters characterising a general unitary $R$.
To avoid searching solutions for such a system, we try a different approach to the problem.
Let us denote with $\tilde{\mathcal H}$ a generator with the required properties:
one that generates the same unitary as $\mathcal H_{\on{CNOT}}$ and contains only 1-qubit interaction terms.
Its general form will be:
\begin{equation}
	\tilde{\mathcal H} = h_0 +
		\sum_{\alpha=1}^3 (h_1^\alpha \sigma_1^\alpha + h_2^\alpha \sigma_2^\alpha ).
\end{equation}
As shown in the main text, for $\HTilde$ to correctly generate CNOT, it must commute with the principal generator $\HH_{\on{CNOT}}$.
Imposing this commutativity removes most of the parameters $h_i^\alpha$, leaving us with the following simplified expression:
\begin{equation}
	\tilde{\HH} = h_0 + h_{1,3} Z_1 + h_{2, 1} X_2.
	\label{eq:cnot_Hprime_reduced}
\end{equation}
The only missing step is now to impose the eigenvalues of $\HH_{\on{CNOT}} - \tilde\HH$ to be integer multiples of $2\pi$.
This gives the following system of equations:
\begin{equation}
\begin{aligned}
	&2\pi \nu_1 = (-\pi + 4 h_0 - 4 h_{1,3} - 4 h_{2,1}) / 4, \\
	&2\pi \nu_2 = (+\pi + 4 h_0 + 4 h_{1,3} - 4 h_{2,1}) / 4, \\
	&2\pi \nu_3 = (+\pi + 4 h_0 - 4 h_{1,3} + 4 h_{2,1}) / 4, \\
	&2\pi \nu_4 = (-\pi + 4 h_0 + 4 h_{1,3} + 4 h_{2,1}) / 4,
\end{aligned}
\end{equation}
with $\nu_i \in \mathbb Z$.
The above system can be seen to have no solution for $h_0, h_{1,3}, h_{2, 1}$,
therefore definitively proving that there is no rotation of the degenerate space,
and integer parameters $\nu_i$, that allow to generate the CNOT gate using only one-qubit interactions.

\subsection{Toffoli gate: derivation through conditions}
\label{subsec:toffoli_main}

We show here the details of how, using Eqs. (1), (2) and (3) in the main text, we can obtain a family of Hamiltonian generators for the Toffoli gate, containing only single- and two-qubit interactions.

The Toffoli gate can be written as
\begin{equation}
	\Toff{} = Z_1^+ + Z_1^-[Z_2^+ + Z_2^- (X_3^+ - X_3^-)],
\end{equation}
where we defined $Z_i^\pm\equiv(1\pm Z_i)/2$, and similarly for $X_i^\pm$.
The principal generator of $\Toff{}$ is $\HToff{} = \pi Z_1^- Z_2^- X_3^-$,
that is, highlighting the three-qubit interaction term,
\begin{equation}
\begin{aligned}
	\HToff =
	(\text{1- and 2-qubit terms}) - \pi/8 \, Z_1 Z_2 X_3.
\end{aligned}
\end{equation}
We start by writing down the general parametrisation of an Hamiltonian containing only one- and two-qubit interactions:
\begin{equation}
	\HTildeToff{} = h_0 \mathds1 +
	\sum h_{i,\alpha} \sigma_i^\alpha +
	\sum J_{i,j}^{\alpha,\beta} \sigma_i^\alpha \sigma_j^\beta.
\end{equation}
Assuming for simplicity that $\HTildeToff{}$ does not contain $Y_i$ components, and imposing the commutativity with $\HToff{}$, we get the following expression:
\begin{equation}
\begin{aligned}
	\HTildeToff{} =& h_0\mathds 1 +
		h_3^x X_3 + h_1^z Z_1 + h_2^z Z_2 + \\
		&J_{13}^{zx} Z_1 X_3 + J_{23}^{zx} Z_2 X_3 + J_{12}^{zz} Z_1 Z_2 + \\
		&(J_{13}^{xx} X_1 + J_{23}^{xx} X_2)(1 + X_3) + \\
		&(J_{12}^{zx} X_2 + J_{13}^{zz} Z_3)(1 + Z_1) + \\
		&(J_{12}^{xz} X_1 + J_{23}^{zz} Z_3)(1 + Z_2).
\end{aligned}
\label{eq:HTildeToff_initial}
\end{equation}
As can be directly verified, the above satisfies $[\HTildeToff{}, \HToff{}] = 0$
while also containing only one- and two-qubit interactions, and no Pauli $Y$ matrices.
We could now directly try and find the set of parameters making the eigenvalues of $\HTildeToff - \HToff$ all integer multiples of $2\pi i$, but the associated calculations are made hard by the many parameters involved.
It turns out however that we can make a number of further assumptions on the form of $\HTildeToff$, and still obtain a viable family of solutions.
One such set of assumptions, that leads to a satisfying family of solutions, is
$J_{12}^{xz} = J_{12}^{zx} = 0$,
$J_{13}^{zx} = J_{23}^{zx} = \pi/8$,
$J_{13}^{xx} = J_{23}^{xx}$,
$J_{23}^{zz} = -J_{13}^{zz}$, and
$h_1^z = h_2^z = -\pi / 8$.
With these assumptions~\cref{eq:HTildeToff_initial} becomes
\begin{equation}
\begin{aligned}
	\HTildeToff{} =
		&h_0 \mathds1 + h_3^x X_3 - \pi/8 (Z_1 + Z_2)(1 - X_3) + \\
		&J_{12}^{zz} Z_1 Z_2 + J_{13}^{xx} (X_1 + X_2)(1 + X_3) + \\
		&J_{13}^{zz} (Z_1 - Z_2)Z_3.
\end{aligned}
\end{equation}
Finally, we impose that the generator is diagonal on the first two qubits, that is, $J_{13}^{xx} = 0$.
Using this simplified expression, $\HH_{\on{Toff}}'{}=\HTildeToff{}-\HToff{}$ becomes
\begin{equation}
\begin{aligned}
	\HH_{\on{Toff}}'{} =
		&\pi/8\, Z_1 Z_2 X_3 + (h_0 - \pi/8)\mathds1 +\\
		&(h_3^x + \pi/8) X_3 + (J_{12}^{zz} - \pi/8) Z_1 Z_2 +\\
		&J_{13}^{zz} (Z_1 - Z_2)Z_3.
\end{aligned}
\label{SM:eq:toffoli_HPrime_reduced}
\end{equation}
With the above simplified expression it is now possible to directly solve the eigenvalue problem.
This results in the following family of solutions:
\begin{equation}
\begin{aligned}

	&\tilde{\mathcal H}_{\on{Toff}} =
	\frac{\pi}{8} \bigg[
	1 + 4 \left(\nu_1 + \nu_2 + 2\nu_3 + \sqrt{(\nu_3-\nu_4)^2} \right) \\
	&- (Z_1 + Z_2)(1 - X_3) + X_3 (-2 - 8\nu_1 + 8\nu_2) \\

	&+ 4 Z_1 Z_2 \left(1/4 + \nu_1 + \nu_2 - 2\nu_3 - \sqrt{(\nu_3-\nu_4)^2} \right) \\
	&+ (Z_2 - Z_1) Z_3 \,\,\sqrt{c(\nu_1, \nu_2, \nu_3, \nu_4)}
	\bigg],
\end{aligned}
\label{eq:toff_tilde_general_solution}
\end{equation}
with
\begin{equation}
\begin{split}
	c(\nu_1, \nu_2, \nu_3, \nu_4) =
		&-(1 + 4\nu_1 - 4\nu_2 + 4\nu_3 - 4\nu_4)\\
		&\times 	(1 + 4\nu_1 - 4\nu_2 - 4\nu_3 + 4\nu_4) =\\
		&\!\!= -[(1+4\nu_{12})^2-(4\nu_{34})^2],
\end{split}
\end{equation}
for all integer values of $\nu_i$ such that $c(\nu_1, \nu_2, \nu_3, \nu_4) \ge 0$.
The corresponding spectrum of $\HH_{\on{Toff}}' = \HTildeToff - \HToff$ is
\begin{equation}
\begin{aligned}
	\lambda_1 &= \lambda_2 = 2\pi \nu_1, \\
	\lambda_3 &= \lambda_4 = 2\pi \nu_2, \\
	\lambda_5 &= \lambda_6 = 2\pi \nu_3, \\
	\lambda_7 &= \lambda_8 = 2\pi (\nu_3 + \lvert\nu_3 - \nu_4\rvert),
\end{aligned}
\end{equation}
while the spectrum of $\HTildeToff$ changes only in that
$\lambda_2 = 2\pi(\nu_1 + 1/2)$.
Consistently with this, $\lambda_2$ is also the eigenvalue corresponding to the non-degenerate eigenspace of $\HToff$, while all the other eigenvalues correspond to eigenvectors orthogonal to this one.
More specifically, we have
\begin{equation}
\begin{aligned}
	\ket{\lambda_1} &= \ket{0, 0, -}, \\
	\ket{\lambda_2} &= \ket{1, 1, -}, \\
	\ket{\lambda_3} &= \ket{1, 1, +}, \\
	\ket{\lambda_4} &= \ket{0, 0, +}, \\
	\ket{\lambda_5} &= \ket{1, 0} \otimes N_5\left[(a - b) \ket0 + \ket1 \right], \\
	\ket{\lambda_6} &= \ket{0, 1} \otimes N_6\left[(a + b) \ket0 + \ket1 \right], \\
	\ket{\lambda_7} &= \ket{1, 0} \otimes N_6\left[(a + b) \ket0 - \ket1 \right], \\
	\ket{\lambda_8} &= \ket{0, 1} \otimes N_5\left[(a - b) \ket0 - \ket1 \right],
\end{aligned}
\label{eq:toffoli_eigenvectors_solution}
\end{equation}
where
\begin{equation}
\newcommand{1 + \bar\nu_{12}}{1 + \bar\nu_{12}}
\begin{gathered}
	a = \frac{\lvert\bar\nu_{34}\rvert}{1 + \bar\nu_{12}{}},
	\qquad
	b = \frac{\sqrt{\bar\nu_{34}^2 - (\bar\nu_{12} + 1)^2}}{1 + \bar\nu_{12}{}},
\end{gathered}
\end{equation}
It is worth noting that the orthogonality of these eigenvectors follows from the easily verified property of the above coefficients: $a^2 - b^2 = 1$.
Furthermore, we note that $c(\nu_1, \nu_2, \nu_3, \nu_4) \ge 0$ cannot be satisfied unless $\nu_3 \neq \nu_4$.
This in turn, looking at \cref{eq:toffoli_eigenvectors_solution},
reveals that all the solutions are made possible by a non-trivial lifting of the degeneracy of the subspaces $\ketbra{0,1}$ and $\ketbra{1,0}$.
Let us now try to understand how and why the derived $\HH_{\on{Toff}}'{}$ works.
Let us use the notation $P_i \equiv \ketbra{\lambda_i}$,
and consider the projector over the last two eigenvectors.
Highlighting the 3-qubit terms, we find
\begin{equation}
\begin{aligned}
	P_7 \simeq - N_6^2 \frac{Z_1 Z_2}{4} \bigg[
		\big((a + b)^2 - 1\big) \frac{Z_3}{2}
		- (a + b) X_3
	\bigg], \\
	P_8 \simeq - N_5^2 \frac{Z_1 Z_2}{4} \bigg[
		\big((a - b)^2 - 1\big) \frac{Z_3}{2}
		- (a - b) X_3
	\bigg].
\end{aligned}
\end{equation}
The term in the Hamiltonian to which these two projectors contribute is
$2\pi \nu_{3,4} (P_7 + P_8)$, with $\nu_{3,4} = \nu_3 + \lvert \nu_3 - \nu_4\rvert$.
A little algebra reveals that the 3-qubit terms in $P_7 + P_8$ are
\begin{equation}
\begin{aligned}
	- \frac{Z_1 Z_2 Z_3}{8} \left[
		N_6^2\big((a+b)^2 - 1\big)
		+ N_5^2 \big((a-b)^2 - 1\big)
	\right] \\
	+ \frac{Z_1 Z_2 X_3}{4} \left[
		N_6^2(a+b) + N_5^2(a-b)
	\right].
\end{aligned}
\end{equation}
Recalling the definitions of $a,b,N_5,N_6$, we see that the coefficient of $Z_1 Z_2 Z_3$ vanishes, and the resulting expression becomes
\begin{equation}
	P_7 + P_8 = (...) +
	Z_1 Z_2 X_3 \frac{1 + 4(\nu_1 - \nu_2)}{16\lvert \nu_3 - \nu_4\rvert}.
\end{equation}
Substitution of the appropriate values of $\nu_i$ shows that the above term can be used to generate the 3-qubit factor $\pi/8 \,\,Z_1 Z_2 X_3$,
\emph{without introducing additional 3-qubit factors}.
In Box~\ref{tcolorbox:toffoli} are given the full expressions for the projectors and the found solutions for the Toffoli gate.
It is also interesting to note that all of the above still holds if the $X_i$ operators are replaced with $Y_i$ operators.
That is, the expressions found solving for the Toffoli, by simple substitution $X_i \to Y_i$,
also give a generator with only 2-qubit interactions for the CCY gate
(that is, the gate that applies $Y$ to the third qubit conditionally to the first 2 qubits being in the $\ket1$ state).

\begin{tbox}[label=tcolorbox:toffoli]{Toffoli}
	\begin{equation*}
	\begin{aligned}
		P_1 &= Z_1^+ Z_2^+ X_3^-,
		\qquad
		P_2 &= Z_1^- Z_2^- X_3^-,
		\qquad
		P_3 &= Z_1^+ Z_2^+ X_3^+,
		\qquad
		P_4 &= Z_1^- Z_2^- X_3^+,
	\end{aligned}
	\end{equation*}
	%
	\begin{equation*}
	\begin{aligned}
		P_5 &= Z_1^- Z_2^+ \frac{1}{2\lvert\bar{\nu}_{34}\rvert}
		\left[
			\lvert\bar\nu_{34}\rvert +
			(1 + \bar\nu_{12}) X_3 -
			\sqrt{-(1 + \bar\nu_{12})^2 + \bar\nu_{34}^2} Z_3
		\right], \\
		P_6 &= Z_1^+ Z_2^- \frac{1}{2\lvert\bar{\nu}_{34}\rvert}
		\left[
			\lvert\bar\nu_{34}\rvert +
			(1 + \bar\nu_{12}) X_3 +
			\sqrt{-(1 + \bar\nu_{12})^2 + \bar\nu_{34}^2} Z_3
		\right], \\
		P_7 &= Z_1^- Z_2^+ \frac{1}{2\lvert\bar{\nu}_{34}\rvert}
		\left[
			\lvert\bar\nu_{34}\rvert -
			(1 + \bar\nu_{12}) X_3 +
			\sqrt{-(1 + \bar\nu_{12})^2 + \bar\nu_{34}^2} Z_3
		\right], \\
		P_8 &= Z_1^+ Z_2^- \frac{1}{2\lvert\bar{\nu}_{34}\rvert}
		\left[
			\lvert\bar\nu_{34}\rvert -
			(1 + \bar\nu_{12}) X_3 -
			\sqrt{-(1 + \bar\nu_{12})^2 + \bar\nu_{34}^2} Z_3
		\right].
	\end{aligned}
	\end{equation*}
	%
	\begin{equation*}
		P_1 + P_2 = \frac{1}{4} (1 + Z_1 Z_2) (1 - X_3),
		\qquad
		P_3 + P_4 = \frac{1}{4} (1 + Z_1 Z_2) (1 + X_3).
	\end{equation*}
	%
	\begin{align*}
		P_5 + P_6 = \frac{1}{4\lvert\bar\nu_{34}\rvert} \left[
			(1 - Z_1 Z_2) \lvert\bar\nu_{34}\rvert +
			(1 - Z_1 Z_2) X_3 (1 + \bar\nu_{12}) +
			(Z_1 - Z_2)Z_3 \sqrt{\bar\nu_{34}^2 - (1 + \bar\nu_{12})^2}
		\right], \\
		P_7 + P_8 = \frac{1}{4\lvert\bar\nu_{34}\rvert} \left[
			(1 - Z_1 Z_2) \lvert\bar\nu_{34}\rvert -
			(1 - Z_1 Z_2) X_3 (1 + \bar\nu_{12}) -
			(Z_1 - Z_2)Z_3 \sqrt{\bar\nu_{34}^2 - (1 + \bar\nu_{12})^2}
		\right].
	\end{align*}

	It is easily verified from the above that
	\begin{equation*}
		P_1 + P_2 + P_3 + P_4 = \frac{1}{2} (1 + Z_1 Z_2),
		\qquad
		P_5 + P_6 + P_7 + P_8 = \frac{1}{2} (1 - Z_1 Z_2),
	\end{equation*}
	so that the sum of the projectors gives the identity as it should.
	On the other hand, multiplying by the appropriate $\nu_i$ factors, we get
	\begin{equation*}
	\begin{aligned}
		2\pi\left[\nu_1(P_1 + P_2) + \nu_2(P_3 + P_4)\right]
		&= (...) + \frac{\pi}{2} (\nu_2 - \nu_1) Z_1 Z_2 X_3, \\
		%
		2\pi\left[\nu_3(P_5 + P_6) + \nu_4(P_7 + P_8)\right]
		&= (...) + \frac{\pi}{2} (\nu_1 - \nu_2) Z_1 Z_2 X_3 + \frac{\pi}{8} Z_1 Z_2 X_3,
	\end{aligned}
	\end{equation*}
	with the last identity holding for $\nu_3\neq \nu_4$.
\end{tbox}

A different way to understand $\HTildeToff$ is to analyse the four two-dimensional subspaces on the main diagonal, exploiting the fact that $\HTildeToff{}$ acts diagonally on the first two qubits.
Straightforward calculations lead to
\begin{equation*}
\begin{aligned}
	\mel{00}{\HTildeToff}{00} &= \pi \left[(\nu_1+\nu_2)-(\nu_1-\nu_2)X\right], \\
	\mel{01}{\HTildeToff}{01} &= 2\pi\nu_3 + \pi\lvert\nu_{34}\rvert(1 - \sigma_{01}), \\
	\mel{10}{\HTildeToff}{10} &= 2\pi\nu_3 + \pi\lvert\nu_{34}\rvert(1 - \sigma_{10}), \\
	\mel{11}{\HTildeToff}{11} &= \frac{\pi}{2}\left[
		(1+2(\nu_1+\nu_2))-(1+2\nu_{12})X
	\right], \\
\end{aligned}
\end{equation*}
where
\begin{equation}
\begin{aligned}
	\sigma_{01}&\equiv\frac{(1+4\nu_{12})X + \sqrt{c}Z}{4\lvert\nu_{34}\rvert}, \\
	\sigma_{10}&\equiv\frac{(1+4\nu_{12})X - \sqrt{c}Z}{4\lvert\nu_{34}\rvert}.
\end{aligned}
\end{equation}
It can be verified that for all values of $\nu_1, \nu_2, \nu_3, \nu_4$, the two-dimensional identity and $X$ are correctly generated in the $\ket{00}$ and $\ket{11}$ spaces, respectively.
On the other hand, in the $\ket{01}$ and $\ket{10}$ spaces, the two-dimensional identity is generated as long as $\nu_3\neq\nu_4$, as was also derived before.

In particular, the class of solutions given by $\nu_1 = \nu_2 = \nu_3 = 0$ is
\begin{equation}
\begin{split}
	\frac{\pi}{8} \bigg[
	1 + 4\lvert\nu_4\rvert
	- 2 X_3 - Z_1 - Z_2
	+ (Z_1 + Z_2) X_3 \\
	+ Z_1 Z_2 (1 - 4\lvert\nu_4 \rvert)
	+ (Z_2 - Z_1) Z_3 \sqrt{16\nu_4^2 - 1}
	\bigg],
\end{split}
\label{SM:eq:toffoli_generator_nu4}
\end{equation}
for all $\nu_4\neq 0$.
It is interesting to look at the explicit form of the exponentials
generated by this class generators.
Computing $\exp(i \tilde{\mathcal H} t)$ with $\tilde{\mathcal H}$ as in \cref{SM:eq:toffoli_generator_nu4}, we get the following unitary:
\begin{equation}
	\begin{pmatrix}
		\mathds1_2 & \mathbb0 & \mathbb0 & \mathbb0 \\
		\mathbb0 & S(t, \nu_4) & \mathbb0 & \mathbb0 \\
		\mathbb0 & \mathbb0 & S(t, \nu_4) & \mathbb0 \\
		\mathbb0 & \mathbb0 & \mathbb0 & X(t)
	\end{pmatrix},
\end{equation}
where
\begin{equation}
	S(t, \nu_4) = \begin{pmatrix}
		a + b & c \\
		c & a - b
	\end{pmatrix},
\end{equation}
\begin{equation}
\begin{gathered}
	a = \frac{1 + e^{2i\pi t \nu_4}}{2},
	\qquad c = \frac{1 - e^{2i\pi t \nu_4}}{8\nu_4},  \\
	b = \frac{(-1 + e^{2i\pi t \nu_4})\sqrt{16\nu_4^2 - 1}}{8\nu_4},
\end{gathered}
\end{equation}
and
\begin{equation}
	X(t) = \frac{1}{2} \begin{pmatrix}
		(1 + e^{i\pi t}) & (1 - e^{i\pi t}) \\
		(1 - e^{i\pi t}) & (1 + e^{i\pi t})
	\end{pmatrix}
\end{equation}
For large (in modulus) values of $\nu_4$,
$a + b \to e^{2i\pi t \nu_4}$, $a - b \to 1$ and $c\to0$,
so that the exponential becomes
\begin{equation}
	\begin{pmatrix}
		\mathds1 & & & & &\\
		& e^{2i\pi t \nu_4} & & & & \\
		& & 1 & & & \\
		& & & e^{2i\pi t \nu_4} & & \\
		& & & & 1 & \\
		& & & & & X(t)
	\end{pmatrix},
\end{equation}
which very closely resembles the matrix obtained by exponentiating the principal generator
$\mathcal H_{\on{Toff}} = \pi Z_1^- Z_2^- X_3^-$:
\begin{equation}
	\exp(i t \mathcal H_{\on{Toff}}) =
	\begin{pmatrix}
		\mathds1 & & & \\
		& \mathds1 & & \\
		& & \mathds1 & \\
		& & & X(t)
	\end{pmatrix}.
\end{equation}
A different solution derived from~\cref{eq:toff_tilde_general_solution} is
\begin{equation}
\begin{aligned}
	\tilde{\mathcal H}_{\on{Toff}} =
	\frac{9\pi}{8} + \frac{3\pi}{4} X_3 - \frac{\pi}{8} (Z_1 + Z_2)
	+ \frac{\pi}{8} Z_1 Z_2 \\
	+ \frac{\pi}{8} (Z_1 + Z_2) X_3
	- \frac{\sqrt7 \pi}{8} (Z_1 - Z_2) Z_3.
\end{aligned}
\end{equation}
Moreover, it is worth noting that~\cref{eq:toff_tilde_general_solution} is only one possible family of solutions, and that different assumptions will lead to different ones.
For example, a similar reasoning as above, starting however from the assumptions $J_{23}^{zz}=J_{13}^{zz}$ will lead to solutions such as (note the use of $(Z_1+Z_2)$ terms here, making this solution not derivable from~\cref{eq:toff_tilde_general_solution}):
\begin{equation}
\begin{aligned}
	\HTildeToff{} =\,
		&\frac{9\pi}{8} - \frac{7\pi}{8}X_3 + \frac{\sqrt{15}\pi}{8} Z_3 + \frac{\pi}{8}Z_1 Z_2 \\
		&\frac{\pi}{8}(Z_1+Z_2)\left(-1+\frac{5}{2}X_3 + \frac{\sqrt{15}}{2}Z_3\right).
\end{aligned}
\end{equation}

\subsection{Toffoli gate: an example of direct \emph{a posteriori} derivation}
\label{subsec:toffoli_posteriori_derivation}
We will here show a line of thought that could have conceivably led to \cref{SM:eq:toffoli_generator_nu4} (in the case $\nu_4=1$), by direct analysis, and without using any of the tools shown in the paper.
It will be useful to keep in mind the expressions of $Z_1 \pm Z_2$:
\begin{equation}
\begin{aligned}
	Z_2 + Z_1 &= \on{diag}(2, 2, 0, 0, 0, 0, -2, -2), \\
	Z_2 - Z_1 &= \on{diag}(0, 0, -2, -2, 2, 2, 0, 0).
\end{aligned}
\end{equation}
Given that we want to generate a CC-X gate, and remembering that
$\exp\left[\frac{i\pi}{2}(1-X)\right] = X$, it is reasonable to start building our Hamiltonian as
\begin{equation}
	\mathcal H_1 = - \pi \left(\frac{Z_1+Z_2}{2}\right)\left(\frac{1-X_3}{2}\right),
\end{equation}
which however will generate an $X$ evolution both in the $\ket{00}$ and in the $\ket{11}$ sectors, while we want it only in the latter sector:
	$\mathcal H_1 \doteq \on{diag}(-X^-, 0, 0, X^-)$.
We can remove the term in the $\ket{00}$ sector exploiting the sign difference introduced by $Z_1 + Z_2$, by directly adding an appropriate 1-qubit interaction term:
\begin{equation}
\begin{aligned}
	\mathcal H_2 &= \frac{1}{2} \left[ \mathcal H_1 + \frac{\pi}{2} (1 - X_3) \right] \\
				 &= \pi \on{diag}(0, X^-\!/2, X^-\!/2, X^-),
\end{aligned}
\label{eq:toffoli_analytical_step2}
\end{equation}
where we remember that $\exp(i\pi X^-) = X$.
\Cref{eq:toffoli_analytical_step2} now correctly reproduces the evolution on $\ket{00}$ and $\ket{11}$, but also wrongly evolves $\ket{01}$ and $\ket{10}$.
To remove these additional terms while at the same time leaving the others unaffected we use the fact that
$\exp(i\pi (1 \pm \bs{\sigma})) = \mathds1$, for any normalised vector of sigma matrices: $\bs\sigma \equiv \sum_{i=1}^3 n_i \sigma_i$ with $n_1^2 + n_2^2 + n_3^2 = 1$.
To convert the central terms in~\cref{eq:toffoli_analytical_step2} into something like this we observe that we can rewrite the second term in the above equation as
\begin{equation}
	\pi/4(1-X_3) = \pi/8(2 - 2 X_3) = \pi/8(5 - 3 - 2 X_3).
\end{equation}
Remembering that $Z_1 Z_2 = \on{diag}(1, -1, -1, 1)$, we substitute the above with
$\pi/8(5 - 3Z_1 Z_2 - 2X_3)$.
This change affects only the central terms, converting the expression into:
$\pi \on{diag}(0, 1 - X/4, 1 - X/4, X^-)$.
The reason this form is preferable is that we can now simply add a factor in the central terms to convert them into an expression of the form $1 - \bs\sigma$.
Adding an interaction of the form $\pi \alpha (Z_2 - Z_1)Z_3$ gives
\begin{equation}
	\pi\on{diag}(0, 1 - X/4 - 2\alpha Z, 1 - X/4 + 2\alpha Z, X^-).
\end{equation}
For the central terms to exponentiate to the identity we need them to become of the form $1-\bs\sigma$ with normalised $\bs\sigma$.
This is easily achieved by choosing $\alpha=\pm\sqrt{15}/8$.
The final expression is therefore:
\begin{equation}
\begin{aligned}
	8/\pi \,\, \mathcal H_3 = -(Z_1 + Z_2)(1 - X_3) \\+ (5 - 3Z_1Z_2 - 2 X_3)
				 \pm \sqrt{15} (Z_2 - Z_1)Z_3.
\end{aligned}
\end{equation}
Note that instead of $\pi \alpha (Z_2 - Z_1)Z_3$ we could have equivalently chosen
$\pi\alpha(Z_2-Z_1)O_3$ for any $O_3 = a Y_3 + b Z_3$ and $a^2 + b^2 = 1$.
The above reasoning explains the origin of the weird $\sqrt{15}$ factor: it comes as the coordinate necessary to make the vector unitary: for $X/4 + x Z / 4$ to be normalized, $x = \sqrt{15}$ must be satisfied.

\section{Supervised learning approach}
\label{sec:supervised_learning_approach}
We here study more in depth the following problem: given a target gate $\G$ and a parametrised Hamiltonian ${\mathcal H(\bs\lambda) = \sum_k \lambda_k \sigma_k}$, with $\lambda_k\in\mathbb R$ and $\sigma_k$ Hermitian operators, what is the set of parameters $\bs\lambda_0$ such that $\exponential(i\HH(\bs\lambda_0)) = \G$?
We present a supervised learning approach to numerically solve this problem, considerably extending the ideas presented in Ref.~\cite{banchi2016quantum}.
Thanks to a number of numerical optimisation techniques, and in particular the use of \ac{AD}~\cite{baydin2015automatic,bartholomewbiggs2000automatic,wengert1964a,bischof2008advances}, we can explore a variety of different scenarios, optimising over potentially hundreds of Hamiltonian parameters.
On top of this, condition (1b) in the main text is used to further speed-up the numerical training, removing many interaction parameters that are known not to lead to the target gate.

\subsection{Supervised learning}
\label{subsec:supervised_learning}

Supervised learning is the task of inferring or approximating a function, given a set of pre-labeled data~\cite{bishop2006pattern,mohri2012foundations}.
A supervised learning algorithm starts with some \emph{model} --- a functional relation $g_{\bs\lambda}$ parametrised by a set of parameters $\bs\lambda$ --- and finds a $\bs\lambda_0$ making $g_{\bs\lambda_0}$ as close as possible to a target function $f$.
To do this, a set of pre-labeled \emph{training data} $\{ (x_1, y_1), (x_2,y_2), ...\}$ is used,
where here $y_k=f(x_k)$ is the output that we want the algorithm to associate to the input $x_k$.

Among the most used supervised learning models are \acp{NN}~\cite{hechtnielsen1989theory,haykin1998neural}.
These are parametric non-linear models which play a prominent role in many machine learning tasks, such as dimensionality reduction, classification, and feature extraction~\cite{hechtnielsen1989theory,haykin1998neural}.
\acp{NN} have also recently proven useful for several problems in quantum many-body theory~\cite{amin2016quantum,wang2016discovering,hush2017machine,carleo2017solving,carrasquilla2017machine,torlai2017manybody,broecker2017quantum,deng2017quantum},
quantum compilation~\cite{swaddle2017generating}, quantum stabilizer codes~\cite{krastanov2017deep} and entanglement quantification~\cite{gray2017measuring}.

A \ac{NN} is \emph{trained} by optimising its parameters using a dataset of pre-labelled data.
A common  way to do this is use variations of a gradient-descent-based technique named \ac{SGD}.
Gradient descent algorithms aim to optimize a problem function $f(\bs x)$, starting from an initial point $\bs x_0$ and performing a number of small steps towards the direction of maximum slope (that is, $\grad f(\bs x)$).
The optimal point ${\bs x}_{\on{opt}}$ is thus obtained via a sequence of small perturbations of the point ${\bs x}$, which starting from ${\bs x}_0$ reaches the nearest local stationary point by following the slope.
In the simplest version of the algorithm, the update rule is simply $\bs x \to \bs x - \eta \grad f(\bs x)$, with $\eta$ a small real parameter commonly referred to as \emph{learning rate}.
\ac{SGD}, on the other hand, is suitable for a situation in which one is given a parametrised functional relationship of the form $f(\bs x; \bs w)$, and asked for a set of ``parameters'' $\bs w_0$ such that $f(\bs x; \bs w_0)$ is minimum (maximum) for all \emph{inputs} $\bs x$.
Such a case can be handled via \ac{SGD}, which in its simplest form involves picking a random $\bs x_1$, executing a number of gradient descent iterations over $\bs w$, then picking a new $\bs x_2$ and iterating the procedure.
The updating rule for \ac{SGD} is therefore of the form
\begin{equation}
	\bs w \to \bs w - \eta\grad_{\bs w}f(\bs x;\bs w).
	\label{eq:updating_rule}
\end{equation}
While standard gradient descent, being a local optimisation algorithm, is liable to getting stuck in local minima, \ac{SGD} can at least partially avoid this issue, in that generally a local minimum for an input $\bs x$ is not a local minimum for a different input $\bs x'$.
Many variations of \ac{SGD} are used in different circumstances.
For example, in the so-called \emph{mini-batch} \ac{SGD}, instead of updating with a single input $\bs x$, one uses a \emph{batch} of inputs $\{\bs x_1,...,\bs x_M\}$, and updates the parameters using the averaged gradient:
$\bs w\to\bs w-\eta\sum_{k=1}^M\grad_{\bs w}f(\bs x_k; \bs w)/M$.
More sophisticated updating rules are used to increase the training efficiency in different circumstances.
Common techniques involve dynamically updating the learning rate, or using \emph{momentum gradient descent}~\cite{goh2017momentum,ruder2016overview} techniques.

To see how this class of optimisation problems is relevant to us, consider the fidelity function $\F$ defined as
\begin{equation}
	\F(\bs\lambda, \psi) \equiv \mel{\psi}{\G^\dagger \exp(i \HH(\bs\lambda)) \G}{\psi},
	\label{eq:def_fidelity}
\end{equation}
with $\G$ the target gate, $\bs\lambda$ the set of parameters, and $\psi$ an input state.
The gate design problem is then equivalent to finding $\bs\lambda$ such that $\F(\bs\lambda, \psi)$ is maximised (that is, equal to 1) for all $\psi$.
One possibility to solve this problem is to consider the average fidelity $\bar{\F}(\bs\lambda)$, for which explicit formulas are known~\cite{banchi2011nonperturbative,pedersen2007fidelity,magesan2011gate}.
Standard optimisation methods, like standard gradient descent or differential evolution, can be applied directly on $\bar{\F}(\bs\lambda)$.
This, however, reveals to be impractical, due to the complexity of the associated parameter landscapes.
On the other hand, \ac{SGD} allows to use a simple and efficient local maximisation method, while at the same time being less prone to getting stuck in local maxima.
This works particularly well in this case, because we know that the sets of parameters corresponding to the target gate are \emph{all and only} those such that \emph{for all inputs $\psi$} the fidelity is unitary.

A crucial step, efficiency-wise, in gradient descent algorithms, is the evaluation of the gradient.
Numerically approximating the gradient, as done in previous works~\cite{banchi2016quantum}, is generally inefficient and scales badly with the number of optimised parameters.
Here we will instead make use of the powerful technique of \acf{AD}~\cite{bartholomewbiggs2000automatic,bischof2008advances}, described in~\cref{subsec:backpropagation}.
\ac{AD} dramatically improves the training efficiency, allowing to explore a richer variety of circumstances.

\subsection{Backpropagation}
\label{subsec:backpropagation}
The gradient evaluation phase is efficiency-wise crucial for the training of a neural network.
Computing the partial derivatives of the cost function with a standard method, like finite differences, has a complexity $\mathcal O(N_{\bs w}^3)$, with $N_{\bs w}$ the number of parameters to differentiate~\cite{bishop2006pattern}.
This inefficiency can however be avoided using \emph{error backpropagation} via \ac{AD}.
With this technique, the complexity of the gradient evaluation phase can be cut down to $\mathcal O(N_{\bs w}^2)$~\cite{bishop2006pattern}.
This works by first decomposing the cost function of the model in terms of elementary operations, that is, functions the gradient of which is known analytically.
In this way the \emph{computational graph} representing the functional relation between input and output is built.
A computational graph is a directed acyclic graph, whose nodes represent the operations, and edges the flowing direction of inputs into outputs (see~\cref{fig:automatic-differentiation}).
Once the computational graph is built, the gradients with respect to the model parameters can be computed efficiently.
This happens in two stages, as schematically illustrated in~\cref{fig:automatic-differentiation}.
At first, every node of the computational graph is progressively computed, starting from the inputs (the current values of the model parameters) up to the final value of the error function.
During this process, the intermediate values of the elementary operations are cached.
This is the so-called \emph{feed-forward} phase.
The second phase (so-called \emph{backpropagation} phase) starts from the output, and consists of progressively computing the gradients of the error function with respect to the independent variables.

To better understand \ac{AD}, let us consider a simple example.
Suppose the error function of the model is of the form
$g(\bs w)\equiv f(\bs f^{(2)}(\bs f^{(1)}(\bs w)))$,
where $\bs w$ is a set of parameters, and $\bs f^{(i)}$ are intermediate ``elementary'' functions, the gradients of which are supposed to be known analytically.
Making use of the chain rule, the gradient of $g$ reads
\begin{equation}
\begin{split}
	\grad g(\bs w) =
	\sum_{k} \partial_{k} f(\bs y^{(2)})
	\grad f_{k}^{(2)}(\bs y^{(1)}),
\end{split}
\label{eq:grad_nested_f}
\end{equation}
where $\bs y^{(2)} = \bs f^{(2)}(\bs f^{(1)}(\bs w))$ and
$\bs y^{(1)} = \bs f^{(1)}(\bs w)$.
During the feed-forward phase the values of $\bs y^{(1)}$ and then $\bs y^{(2)}$ are progressively computed and cached.
Using $\bs y^{(2)}$, and the known expression for $\partial_k f$, $\partial_{k} f(\bs y^{(2)})$ is then efficiently computed.
The process continues by evaluating $\grad f_k^{(2)}$, which is written as
\begin{equation}
	\grad f_k^{(2)}(\bs y^{(1)}) =
	\sum_j \partial_j f_k^{(2)}(\bs y^{(1)}) \grad f_j^{(1)}(\bs w).
\end{equation}
Again, being $\bs y^{(1)}$ already computed during the feed-forward, $\partial_j f_k^{(2)}(\bs y^{(1)})$ is readily computed.
The last component needed for the full gradient is $\partial_i f_j^{(1)}(\bs w)$, all parts of which are known.
This method therefore allows to efficiently evaluate numerical the gradient of complicated functions, without approximating the derivatives.

In the context of training neural networks, the function to be derived is the \emph{cost function} of the network, that is, roughly speaking, the (euclidean) distance between the result obtained for an input and the corresponding training output.
For the gate design problem, we will use another notion of \emph{distance} between output obtained and output expected.
For quantum states, the fidelity between these turns out to work well.

\begin{figure*}[tb]
	\centering
	\begin{minipage}{0.3\linewidth}
		\centering
		\includegraphics[height=1.3\textwidth]{automatic-differentiation}
	\end{minipage}\hfill
	\begin{minipage}{0.3\linewidth}
		\centering
		\includegraphics[height=1.3\textwidth]{automatic-differentiation2}
	\end{minipage}\hfill
	\begin{minipage}{0.3\linewidth}
		\centering
		\includegraphics[height=1.3\textwidth]{automatic-differentiation3}
	\end{minipage}
	\caption{
		Examples of \ac{AD} in backpropagation mode.
		\textbf{(a)}
		Schematic representation of \ac{AD} of a function with one output and two inputs.
		Starting from numerical values for $x_1$ and $x_2$, one computes $g(x_1, x_2)$ and then $f(g(x_1, x_2))$.
		To get $\grad f(g(x_1, x_2))$, one then computes $f'(g(x_1, x_2))\partial_i g(x_1, x_2)$.
		Note that all components of this expression are known: $f'$ and $\partial_i g$ are known by assumption, and the value of $g(x_1, x_2)$ has been computed and cached in the forward propagation phase.
		\textbf{(b)} Example of application of \ac{AD} to compute the gradient of $\cos(x_1 x_2)$.
		\textbf{(c)} Using the same example function as (b), we give an example of the actual number computed at all stages when the inputs are $(x_1, x_2) = (\pi / 2, 2)$.
	}
	\label{fig:automatic-differentiation}
\end{figure*}

\subsection{Implementation details}
\label{subsec:implementation_details}

We used python as language of choice for the implementation of the supervised learning.
Being Python language of widespread use in the machine learning community, many libraries and frameworks are available to build computational graphs over which \ac{AD} can be used.
In particular, we used \textsc{Theano}~\cite{team2016theano}, together with the \textsc{QuTiP} library for the simulation of the dynamics of quantum systems~\cite{johansson2012qutip,johansson2013qutip}.

Our implementation allows the training of an arbitrary target gate, parametrised via a time-independent Hamiltonian $\mathcal H(\bs\lambda)$.
The parametrisation is completely arbitrary (provided the dependence on the parameters is linear), so that the Hamiltonian can be chosen as $\mathcal H(\bs\lambda)=\sum_i \lambda_i A_i$ for any set of matrices $A_i$ and number of parameters $\lambda_i$.
This is made possible by the flexibility of \ac{AD}, which allows to automatically build an efficiently differentiable computational graph, without needing to hardcode the structure of the Hamiltonian.

The goal of the algorithm is, given a target gate $\mathcal G$ and a parametrisation for the Hamiltonian $\mathcal H(\bs\lambda)$, find the $\bs\lambda_0$ such that $\exp(i\mathcal H(\bs\lambda_0))=\mathcal G$.
We use for the purpose mini-batch \ac{SGD} with momentum.
The \emph{mini-batch} version of \ac{SGD} involves computing the gradient, at every iteration, averaging over the gradients computed for a number of states.
Making such \emph{batches} of states larger or smaller allows to enhance or decrease the variance of the gradients with respect to the input state.
The use of \emph{momentum}~\cite{ruder2016overview,goh2017momentum} involves using a modified version of~\cref{eq:updating_rule}.
The updating rule is instead given by
\begin{equation}
\begin{aligned}
	\bs v &\to \gamma \bs v + \eta \grad_{\bs\lambda} \F(\psi, \bs\lambda), \\
	\bs\lambda &\to \bs\lambda + \bs v,
\end{aligned}
\label{eq:updating_rule_momentum}
\end{equation}
where here $\eta$ is the \emph{learning rate} and $\gamma$ the \emph{momentum}.
The use of the auxiliary parameter $\bs v$ during the training discourages sudden changes of direction, and can make the training significantly more efficient~\cite{goh2017momentum}.

While the cost function $\F$ is always real, some of the intermediate calculations needed to compute it involve complex numbers.
While this poses no fundamental problems, many of the \ac{ML} libraries do not support \ac{AD} over functions with complex inputs or outputs.
We worked around this problem using a similar trick to the one reported in~\cite{leung2017speedup}.
In particular, to use the existing framework, we mapped the problem into one involving only real numbers.
To do this, we map complex matrices into real ones via the bijection
$A\mapsto\mathfrak{Re}(A)\equiv\mathds1\otimes A_{R} - i \sigma_y\otimes A_{I}$,
where $A_R$ and $A_I$ are the real and imaginary parts of $A$, respectively.
At the same time, state vectors are to be mapped to
$\Psi\mapsto\mathfrak{Re}(\Psi)\equiv(\Psi_R, \Psi_I)^T$.
It is easy to verify that with this mapping
$A\Psi\mapsto\mathfrak{Re}(A\Psi)=\mathfrak{Re}(A)\mathfrak{Re}(\Psi)$,
so that all calculations can be equivalently be carried out with the real versions of matrices and vectors.

More specifically, the employed algorithm involves the following steps:
\begin{enumerate}
	\item Choose an initial set of parameters $\bs\lambda$ (randomly, or specific values if one has an idea of where a solution might be).
	A number of other hyperparameters have to be decided at this step, depending on the exact \ac{SGD} method used. In particular, for mini-batch \ac{SGD} with momentum and decreasing learning rate, one has to decide the momentum $\gamma$, the initial value of $\eta$, the rate at which $\eta$ decreases during the training, and the size $N_b$ of the batches of states used for every gradient descent step.
	\item Repeat the following loop $N_e$ times, or until a satisfactory result is obtained.
	Each such iteration is conventionally named an \emph{epoch}.
	Another hyperparameter to be chosen beforehand is the number of training states $N_{tr}$ to be used in each epoch.
	Once this is fixed, every epoch will involve a number $N_{tr}/N_e$ of gradient descent steps, each one using $N_e$ states for a single gradient calculation.
	$N_e$ random training states are sampled, to be used during the epoch.
	\begin{enumerate}
		\item Pick $N_b$ of the $N_e$ training states.
		\item Forward-propagate each state of the sample, and then backpropagate the gradients, thus computing the average gradient over the mini-batch $\grad_{\bs\lambda} \F(\bs\lambda)$.
		\item Update the coupling strengths $\lambda$ as per~\cref{eq:updating_rule_momentum}.
		\item Return to point (a).
	\end{enumerate}
\end{enumerate}

\subsection{Results}
\label{subsec:numerical_results}

A sample of training results for Toffoli, Fredkin, and "double Fredkin" gates are given in Fig. 1 (a), (b), and (c) in the main text.
In~\cref{fig:toffoli_diagonal_parhistories,fig:fredkin_diagonal_parhistories,fig:doublefredkin_diagonal_parhistories} are shown the training histories of the parameters for eight different solutions for Toffoli, Fredkin and \emph{double Fredkin}, respectively.
These illustrate how quickly the networks converge for different initial values of the parameters.
In all the shown cases the target gates are obtained with unit fidelity up to numerical precision (that is, all fidelities are between $1-10^{-16}$ and $1$).
Different sets of optimisation hyperparameters are found to give acceptable solutions.
For the trainings shown in this paper we used a dynamically updated learning rate given, for the $k^{\text{th}}$ epoch, by $\eta=1/(1 + k \alpha)$ with the \emph{decay rate} $\alpha=0.005$.
The other hyperparameters were chosen as
$\gamma=0.5$,
$N_b = 2$, $N_{tr} = 200$.
Different initial values for the parameters were tested, but in most cases we started the training with either vanishing or random (following a normal distribution) parameters.
For the training of the four-qubit gate we found the network to converge sooner to a solution when the parameters were initialised to a positive value (often with all parameters initialised to $4$).

In~\cref{fig:toffoli_fidVSpars,fig:fredkin_fidVSpars,fig:doublefredkin_fidVSpars} we report the behaviour of the fidelity upon changes of the learnt Hamiltonian parameters, for Toffoli, Fredkin and \emph{double Fredkin} gates, respectively.
As shown in these plots, the stability of the implemented gates with respect to variations of time and interactions values greatly varies between different solutions, as well as between different parameters in the same solutions.

To assess the feasibility of nontrivial gates in more restrictive experimental scenarios, we performed a systematic analysis of the reachability of Fredkin and Toffoli gates when allowing only for single-qubit and $X_i X_j+Y_i Y_j$ two-qubit interactions, and in the less restrictive setting of allowing for all $X_i X_j$ and $Y_i Y_j$ interactions.
The results are shown in~\cref{fig:fredkin_XY,fig:fredkin_XX,fig:toffoli_XX,fig:toffoli_XY}.
For the Fredkin gate, in the more restrictive $XX$ interactions setting, the biggest fidelity obtained was $\F\simeq 0.94$, while when allowing for all $XX$ and $YY$ interactions the maximum fidelity obtained was $\F\simeq0.999$.
For the Toffoli gate, the maximum fidelity obtained in the $XX$ scenario was $\F\simeq0.94$ as well, while when allowing for all $XX$ and $YY$ interactions the best training results corresponded to $\F\simeq 0.98$.
To have more consistent results, in all the training instances shown here all the hyperparameters, except for the interaction parameters' initial values, were chosen to have the same value.
In particular, each training instance was run for $200$ epochs, each one using $200$ random quantum states as inputs, divided in batches of $2$ elements.
This choice of hyperparameters is mostly empirical, and it is possible for different values to provide better results.

The above provides further evidence for the flexibility of the supervised learning approach, which can produce solutions with good fidelities even in more restrictive scenarios, closer to the capabilities of state of the art experimental architectures.
Furthermore, the values of the interaction strengths for many of the presented solutions are found to be compatible with the capabilities of state of the art circuit-QED architectures with gate times of the order of tens of nanoseconds~\cite{potocnik2018studying}.

Additional solutions and data, as well as the code used to produce them, is available in the GitHub repository
\href{https://github.com/lucainnocenti/quantum-gate-learning-1803.07119}{lucainnocenti/quantum-gate-learning-1803.07119}.
This repository contains all the code used to reproduce the solutions presented in this paper, as well as to train arbitrary gates on arbitrary numbers of qubits.
Even more generally, arbitrary (linearly) parametrised matrices can be used as training model, allowing a high degree of flexibility.

\begin{figure*}[]
	\centering
	\begin{tikzpicture}
		\node[anchor=south west] (A) at (0, 0)%
			{\includegraphics[width=.8\linewidth]{toffoli_diagonal_parhistories}};
		\node[above] at (4., -.2) {$t_e$};
		\node[above] at (11.5, -.2) {$t_e$};
	\end{tikzpicture}%

	\caption{
		Training histories for the \textbf{Toffoli} gate with only diagonal interactions.
		In each plot are reported the values of the $9$ network parameters during the training process, for each training epoch $t_e$.
		Each training process was left running until convergence to unit fidelity, therefore, the number of epochs in the horizontal axes differs for different trainings instances.
		The histories shown here correspond to training instances in which the parameters were initialised at various values, as seen from the leftmost values in each plot.
	}
	\label{fig:toffoli_diagonal_parhistories}
\end{figure*}

\begin{figure*}[]
	\centering
		\begin{tikzpicture}
		\node[anchor=south west] (A) at (0, 0)%
			{\includegraphics[width=.8\linewidth]{fredkin_diagonal_parhistories}};
		\node[above] at (4., -.2) {$t_e$};
		\node[above] at (11.5, -.2) {$t_e$};
	\end{tikzpicture}%

	\caption{
		Training histories for the \textbf{Fredkin} gate with only diagonal interactions.
		In each plot are reported the values of the $9$ network parameters during the training process, for each training epoch $t_e$.
		Each training process was left running until convergence to unit fidelity, therefore, the number of epochs in the horizontal axes differs for different trainings instances.
		The histories shown here correspond to training instances in which the parameters were initialised at various values, as seen from the leftmost values in each plot.
	}
	\label{fig:fredkin_diagonal_parhistories}
\end{figure*}

\begin{figure*}[]
	\centering
	\begin{tikzpicture}
		\node[anchor=south west] (A) at (0, 0)%
			{\includegraphics[width=.8\linewidth]{doublefredkin_diagonal_initvalues4_parhistories}};
		\node[above] at (4., -.2) {$t_e$};
		\node[above] at (11.5, -.2) {$t_e$};
	\end{tikzpicture}%

	\caption{
		Training histories for the \textbf{double-Fredkin} gate with only diagonal interactions.
		In each plot are reported the values of the $18$ network parameters during the training process, for each training epoch $t_e$.
		Each training process was left running until convergence to unit fidelity, therefore, the number of epochs in the horizontal axes differs for different trainings instances.
		All the histories shown here correspond to training instances in which the parameters were initialised to $4$.
	}
	\label{fig:doublefredkin_diagonal_parhistories}
\end{figure*}

\newcommand{\graphicswithlabel}[3]{%
	\begin{tikzpicture}
		\node[anchor=north west] (A) at (0, 0)%
			{\includegraphics[width=.48\linewidth]{#2}};
		\node[above] at (1.4, -.85) {\small\textbf{(#1)}};
		\node[above] at (4.6, -4.65) {#3};
		\node[above] at (1., -.2) {$\F_{\bs{\lambda}}(\psi)$};
	\end{tikzpicture}%
}
\begin{figure*}
	\centering
	\graphicswithlabel{a}{toffoli_fidVStime}{$\alpha$}\vspace{-6pt}
	\graphicswithlabel{b}{toffoli_fidVStime_wide}{$\alpha$}
	\graphicswithlabel{c}{toffoli_fidVSpar0}{$\lambda_1$}\vspace{-6pt}
	\graphicswithlabel{d}{toffoli_fidVSpar1}{$\lambda_2$}
	\caption{
		Fidelity $\F_{\bs{\lambda}}(\psi)$ vs variations of ${\bs{\lambda}}$, for different test states, for the \textbf{Toffoli} gate. The five test states $\psi$ are sampled randomly.
		\textbf{(a)} Global relative variations of ${\bs{\lambda}}$, that is, plotting the fidelity against $\alpha{\bs{\lambda}}$ for $0.9 \le\alpha\le 1.1$.
		Note that this is equivalent to studying how the fidelity changes with respect to uncertainties in the evolution time, that is, how much does $\exp(i\HH t')$ differ from $\exp(i\HH t)$.
		\textbf{(b)} Same as \emph{(a)} but with $0\le\alpha\le 1.2$.
		\textbf{(c)} Plot of $\F_{\bs{\lambda}}(\psi)$ against \emph{absolute} variations of a single element of ${\bs{\lambda}}$, in this case the first one, i.e. we take $\lambda_1\in[-10,10]$.
		\textbf{(d)} Like \emph{(c)} but for $\lambda_2$.
	}
	\label{fig:toffoli_fidVSpars}
\end{figure*}

\renewcommand{\graphicswithlabel}[3]{%
	\begin{tikzpicture}
		\node[anchor=north west] (A) at (0, 0)%
			{\includegraphics[width=.48\linewidth]{#2}};
		\node[above] at (1.4, -.85) {\small\textbf{(#1)}};
		\node[above] at (4.6, -4.65) {#3};
		\node[above] at (1., -.2) {$\F_{\bs{\lambda}}(\psi)$};
	\end{tikzpicture}%
}
\begin{figure*}
	\centering
	\graphicswithlabel{a}{fredkin_fidVStime}{$\alpha$}\vspace{-6pt}
	\graphicswithlabel{b}{fredkin_fidVStime_wide}{$\alpha$}
	\graphicswithlabel{c}{fredkin_fidVSpar0}{$\lambda_1$}\vspace{-6pt}
	\graphicswithlabel{d}{fredkin_fidVSpar3}{$\lambda_2$}
	\caption{
		Fidelity $\F_{\bs{\lambda}}(\psi)$ vs variations of ${\bs{\lambda}}$, for different test states, for the \textbf{Fredkin} gate. The five test states $\psi$ are sampled randomly.
		\textbf{(a)} Global relative variations of ${\bs{\lambda}}$, that is, plotting the fidelity against $\alpha{\bs{\lambda}}$ for $0.9 \le\alpha\le 1.1$.
		Note that this is equivalent to studying how the fidelity changes with respect to uncertainties in the evolution time, that is, how much does $\exp(i\HH t')$ differ from $\exp(i\HH t)$.
		\textbf{(b)} Same as \emph{(a)} but with $0\le\alpha\le 1.2$.
		\textbf{(c)} Plot of $\F_{\bs{\lambda}}(\psi)$ against \emph{absolute} variations of a single element of ${\bs{\lambda}}$, in this case the first one, i.e. we take $\lambda_1\in[-10,10]$.
		\textbf{(d)} Like \emph{(c)} but for $\lambda_2$.
	}
	\label{fig:fredkin_fidVSpars}
\end{figure*}

\renewcommand{\graphicswithlabel}[3]{%
	\begin{tikzpicture}
		\node[anchor=north west] (A) at (0, 0)%
			{\includegraphics[width=.48\linewidth]{#2}};
		\node[above] at (1.4, -.85) {\small\textbf{(#1)}};
		\node[above] at (4.6, -4.65) {#3};
		\node[above] at (1., -.2) {$\F_{\bs{\lambda}}(\psi)$};
	\end{tikzpicture}%
}
\begin{figure*}
	\centering
	\graphicswithlabel{a}{doublefredkin_fidVStime}{$\alpha$}\vspace{-6pt}
	\graphicswithlabel{b}{doublefredkin_fidVStime_wide}{$\alpha$}
	\graphicswithlabel{c}{doublefredkin_fidVSpar0}{$\lambda_1$}\vspace{-6pt}
	\graphicswithlabel{d}{doublefredkin_fidVSpar1}{$\lambda_2$}
	\caption{
		Fidelity $\F_{\bs{\lambda}}(\psi)$ vs variations of ${\bs{\lambda}}$, for different test states, for the \textbf{double Fredkin} gate. The five test states $\psi$ are sampled randomly.
		\textbf{(a)} Global relative variations of ${\bs{\lambda}}$, that is, plotting the fidelity against $\alpha{\bs{\lambda}}$ for $0.9 \le\alpha\le 1.1$.
		Note that this is equivalent to studying how the fidelity changes with respect to uncertainties in the evolution time, that is, how much does $\exp(i\HH t')$ differ from $\exp(i\HH t)$.
		\textbf{(b)} Same as \emph{(a)} but with $0\le\alpha\le 1.2$.
		\textbf{(c)} Plot of $\F_{\bs{\lambda}}(\psi)$ against \emph{absolute} variations of a single element of ${\bs{\lambda}}$, in this case the first one, i.e. we take $\lambda_1\in[-10,10]$.
		\textbf{(d)} Like \emph{(c)} but for $\lambda_2$.
	}
	\label{fig:doublefredkin_fidVSpars}
\end{figure*}

\newcommand{\traininggraphicswithlabel}[1]{%
	\begin{tikzpicture}
		\node[anchor=north west] (A) at (0, 0)%
			{\includegraphics[width=.85\linewidth]{#1}};
		\node[above] at (0, -.8) {$\barF$};
	\end{tikzpicture}%
}

\begin{figure*}
	\centering

	\traininggraphicswithlabel{figures/fredkin_noancillae_XXall.pdf}
	\caption{
		Training results for the Fredkin gate with a generator containing one-qubit interactions and two-qubit interactions of the form $J_{ij}(X_i X_j + Y_i Y_j)$.
		Every point shows the final fidelity obtained at the end of a training procedure.
		The hyperparameters, as well as the total number of training iterations, are kept the same in all the training instances shown here.
		The initial parameters' values are the same within each horizontal sector, but changed between different sectors.
		The initial parameters' values within each sector have been chosen as all equal to $c$ (that is, $\lambda_i=c$ for all $i$). The values of $c$ are $0, 1,..., 10$, with $c=0$ in the leftmost sector and $c=10$ in the last to rightmost one.
		The rightmost sector contains the results of training attempts with the initial values chosen at random (that is, with $\lambda_i$ sampled according to the uniform normal distribution, independently for each $i$).
		The greatest reached fidelity, obtained with initial parameters' values $\lambda_i=3$, is $\F\simeq0.94$.
	}
	\label{fig:fredkin_XX}
\end{figure*}

\begin{figure*}
	\centering
	\traininggraphicswithlabel{figures/fredkin_noancillae_XYall.pdf}
	\caption{
		Training results for the Fredkin gate with a generator containing all one-qubit interactions and two-qubit interactions of the form $J^{(1)}_{ij}X_i X_j + J^{(2)}_{ij}Y_i Y_j$.
		The initial conditions are chosen as in~\cref{fig:fredkin_XX}.
		The maximum fidelities obtained are $\F\simeq0.999$, obtained in multiple instances $c=3$ and $c=4$ sectors.
	}
	\label{fig:fredkin_XY}
\end{figure*}

\begin{figure*}
	\centering
	\traininggraphicswithlabel{figures/toffoli_noancillae_XXall.pdf}
	\caption{
		Training results for the Toffoli gate with a generator containing all one-qubit interactions and two-qubit interactions of the form $J_{ij}(X_i X_j + Y_i Y_j)$.
		The initial conditions are chosen as in~\cref{fig:fredkin_XX}.
		The maximum fidelities obtained are $\F\simeq0.94$, obtained in the $c=3$ sector.
	}
	\label{fig:toffoli_XX}
\end{figure*}

\begin{figure*}
	\centering
	\traininggraphicswithlabel{figures/toffoli_noancillae_XYall.pdf}
	\caption{
		Training results for the Toffoli gate with a generator containing all one-qubit interactions and two-qubit interactions of the form $J^{(1)}_{ij}X_i X_j + J^{(2)}_{ij}Y_i Y_j$.
		The initial conditions are chosen as in~\cref{fig:fredkin_XX}.
		The maximum fidelities obtained are $\F\simeq0.98$, obtained in the $c=4$ and $c=5$ sectors.
	}
	\label{fig:toffoli_XY}
\end{figure*}

\end{document}

redpajama_set_name: RedPajamaArXiv