\section{Introduction}

GRO\,J2058+40 was discovered as a 198 s X-ray pulsar by BATSE on board the
Compton Gamma-Ray Observatory ({\it CGRO}) during a giant outburst in 1995
September-October \citep{wilson98}. This giant (type II) outburst was followed
by a series of weaker (type I) outbursts, whose intensity peaked at intervals of
about 55 days \citep{wilson98,wilson05}. Because BATSE detected the odd
outbursts, counting from the giant outburst, to be brighter than even outbursts,
it was proposed that \gro\ was undergoing periastron and apastron outburst in a
110 day orbit \citep{wilson98}. However, the all-sky monitor (ASM)
on board the {\it Rossi X-Ray Timing Explorer (RXTE)} did not see such a pattern.
Consequently, the orbital period was established to be 55 days \citep{wilson05}.

\gro\ was X-ray active from 1996 to 2002, but only detected at periastron
passages. Pulsations were not detected during the epoch 2003--2004
\citep{wilson05}.
\gro\ entered a new period of enhanced X-ray activity in May 2008
\citep{krimm08}. Although the X-ray flux doubled in less than a week, the
luminosity was relatively low, $L_x=1.5\times 10^{36}$ \ergs\ assuming a
distance of 9 kpc \citep{wilson05,reig15}. This is more consistent with  the lower-luminosity type I outbursts. However, because the 2008 X-ray data were
limited to just one {\it Swift}/XRT pointed observation, it is likely that
the peak of the outburst was missed.

Another giant X-ray outburst occurred in March 2019 with a peak luminosity in
the energy range 3--78 keV of $L_x=5.6 \times 10^{37}$ \ergs\ \citep{kabiraj20},
similar to the discovery luminosity back in 1995.  {\it NuSTAR} and {\it
AstroSat} observations performed during this outburst revealed the possible
presence of a cyclotron line at 10 keV together with some harmonics
\citep{molkov19,mukerjee20}. During {\it AstroSat} observations, a
quasi-periodic oscillation at 0.090 Hz was detected during the decay of the
outburst \citep{mukerjee20}. The {\it NuSTAR} observations only detected the
cyclotron line in a very narrow range of the spin phases of \gro\ and covering
only $\sim10$\% of the entire spin period \citep{molkov19}, but it did not seem
to be present in the pulse-average spectrum \citep{kabiraj20}.

The optical counterpart to GRO\,J2058+40 is an O9.5-B0e IV-V star
\citep{reig05a,wilson05}. The massive companion is a long-term variable system
with photometric and spectroscopic changes on timescales of years
\citep{kiziloglu07b,reig15,reig16}. The underlying B star also displays fast
optical multiperiodic variability that is  attributed to nonradial pulsations
\citep{kiziloglu07b,reig22}. We present the most complete and
detailed study of the optical variability of the BeXB \gro\ performed so far.

\section{Observations and data analysis}

We obtained optical photometry, spectroscopy, and polarimetry of the optical
counterpart to \gro. The observations were made from the 1.3 m telescope at the
Skinakas observatory\footnote{https://skinakas.physics.uoc.gr/en/} (SKO) in
Crete (Greece). Additional spectra were retrieved from the data archive of the
2.0 m Liverpool Telescope\footnote{https://telescope.livjm.ac.uk/}(LT).

\subsection{Photometry}

Photometric observations with the Johnson-Cousins $B$,
$V$, $R,$ and $I$ filters were made with the 1.3 m telescope of the
Skinakas Observatory. The telescope
was equipped with a  2048$\times$2048 ANDOR CCD with a 13.5 $\mu$m pixel
size. In this configuration, the  plate scale is 0.28$\arcsec$/pixel,
hence providing a field of view of $9.5 \times 9.5$ arcmin$^2$. Standard
stars from the Landolt list \citep{landolt09} were used for the
transformation equations.  The data were reduced in the
standard way using the IRAF tools for aperture photometry\footnote{A User's Guide to CCD
Reductions with IRAF, Philip Massey, February 1997. https://iraf.net/irafdocs}
\citep{tody86}. After the
standardization process, we finally assigned an error to the calibrated
magnitudes of the target given by the rms of the residuals between the
cataloged and calculated magnitudes of the standard stars. The photometric
magnitudes are given in Table~\ref{photres} (see Appendix~\ref{app}).

\subsection{Polarimetry}

Polarimetric observations were made with the RoboPol photopolarimeter attached
to the focus of the 1.3 m telescope of the Skinakas Observatory.  In the
polarimetry configuration a plate scale of  0.43 $\arcsec$/pixel is achieved
with a 2048$\times$2048 ANDOR CCD with a 13.5 $\mu$m pixel size.  RoboPol is an
imaging photopolarimeter that simultaneously measures the Stokes parameters of
linear polarization  of all sources in the field of view
\citep{king14,ramaprakash19}. RoboPol splits the incident light in two beams,
each half incident on a half-wave retarder followed by a Wollaston prism.  The
fast axis of the half-wave retarder in front of the first prism is rotated by
$67.5^{\circ}$ with respect to the other retarder. Every point in the sky is
thereby projected to four points on the CCD with different polarization state.
The photon counts in each spot, measured using aperture photometry, were used to
calculate the $U$ and $Q$ parameters of linear polarization.  To optimize the
instrument sensitivity, a mask was placed in the telescope focal plane. The
absence of moving parts allows RoboPol to compute all the Stokes parameters of
linear polarization in one shot. The intrumental polarization and polarization
angle zeropoint were controlled with regular measurements of polarimetric
standards as described in \citet{blinov21}. All observations were corrected for the
instrumental polarization, and uncertainties were propagated accordingly. The
results of the polarimetric observations are given in Table~\ref{polres}
(Appendix~\ref{app}).

\subsection{Spectroscopy}

The 1.3\,m telescope at SKO  was equipped with a
2048$\times$2048 (13.5 micron) pixel ANDOR IKON CCD and a 1302 l~mm$^{-1}$
grating, giving a nominal dispersion of $\sim$0.9 \AA/pixel.  Spectra of
comparison lamps were taken before and after each exposure in order to perform
the wavelength calibration and account for small variations produced by the
tension and flexture of the telescope at different zenith angles during the
night that could affect the wavelength calibration The 2.0 m LT is a fully
robotic telescope at the Roque de Los Muchachos on the Canary Island of La Palma
(Spain). We downloaded and analyzed all the spectra available at the LT data
archive for this source.  These spectra were obtained  with the Fibre-fed
RObotic Dual-beam Optical Spectrograph (FRODOSpec), which is  a dual-beam design
multipurpose integral-field input spectrograph that splits the beam before the
entrance to the individually optimized collimators. We used the
spectra taken with the red arm, which cover the wavelength range 5900--8000 \AA\
with a dispersion of 0.6 \AA/pixel. A Xenon arc exposure was obtained after the
science exposure. The log of the spectroscopic observations is given in
Table~\ref{specres} (see Appendix~\ref{app}).

The \ha\ line is the prime indicator of the state of the disk. The structure and
size of the disk affect not only the strength of the \ha\ line, but also its
shape. We investigated the evolution of the strength and line
profile over about an 18 year period. We took the equivalent width of the \ha\
line (\ew) as a proxy of the strength of the line. To study the changes
experienced by the shape of the line, we fit the line profile with one or two
Lorentzian functions, depending on whether the line displayed a single-peak or a
double-peak profile. In the case of a split profile, the peak at shorter
wavelength is called the blue or violet peak, and it is denoted by $V$, and the
peak at the longer wavelength is called the red peak and is represented by $R$.

The three fitting parameters are the central wavelength ($\lambda_i$), the
full width at half maximum ($FWHM_i$), and the relative intensity of the peak
($I_i$). Here the subindex $i$ refers to either the blue ($V$) or red peak
($R$). When the line showed a double-peak profile, we obtained the
following quantities: $\Delta V=\lambda_{\rm R}-\lambda_{\rm V}$
and $V/R=\log(I_{\rm V}/I_{\rm R})$, where $I_{\rm V}$ and $I_{\rm R}$ are the
intensity of the blue and red peaks relative to the continuum. For single-peak
profiles, $\Delta V=0$ and $\log(V/R)$ was not defined.

The \ew\ was computed by means of a numerical integration method, the
trapezoidal rule, as implemented in the scipy python package. To ensure an
homogeneous processing, all the spectra were normalized with respect to the
local continuum, which  was rectified to unity by employing a polynomial fit.
Because the definition of the continuum is the main source of uncertainty in the
computation of the equivalent width, we averaged over different selections of the
continuum and different polynomial fits (different grade). A total of 24
measurements were obtained for each spectrum. The error in \ew\ is the standard
deviation of these measurements.

For the line fits, we used the {\tt deblend} task of the SPLAT package
\citep{skoda14} of the Starklink project \citep{currie14}. The fits require the
continuum to be at zero level. Therefore, the best polynomial fit of the
continuum was subtracted from the spectra prior to the Lorentzian fits.

\begin{figure}
\resizebox{\hsize}{!}{\includegraphics{./fig1.eps}}
\caption[]{Evolution of the optical observables. From top to bottom, we show the \ha\
equivalent width, continuum emission in the $R$ band, $B-V$ color index, and
polarization degree in the $R$ band (not corrected for ISM polarization).}
 \label{opt}
\end{figure}

\section{Results}

\subsection{Brightness variability}
\label{res}

Figure~\ref{opt} shows the long-term optical variability of \gro. In general,
the correlation among all optical observables is good. As the strength of
the line emission increases, the continuum emission and the polarization
increase.
\gro\ went through two high and three low optical states during the course of
our observations.

A closer look at Fig.~\ref{opt} reveals some interesting results. First, the
duration of the line emission in the low-optical states is significantly shorter
than the continuum emission, that is, the changes in \ha\ occur on faster timescales than in the continuum. While the continuum emission may remain stable for
years during the optical minimum, the \ha\ line emission reaches a minimum and
begins to increase after just a few weeks. Second, the continuum emission reaches
its peak by about $\sim400$ days before the \ew\ does. Third, the source
exhibited stronger line and continuum emission  during the 2019 high state than
during the 2009 high state: \ew$=-13.6$ \AA\ and $R=13.74$ mag compared to
\ew$=-10$ \AA\ and $R=13.92$ mag.

The transition from low to high states and then back to low states occurred
smoothly over a period of several years. Overall, the rate of change in the \ew\ is
roughly similar during the two transitions from the low to the high states as
it is for the two declining phases. However, the changes in the optical
observables during the declining phases are in general
more abrupt that during the rising phases. To quantify this result, we performed a
linear fit to the \ew\ for the four different phases, as shown in
Fig.~\ref{ewfit}. The results are summarized in Table~\ref{linfit}. The slopes of
the fits are 1.5--2 times larger during the decay phase.

\begin{figure}
\resizebox{\hsize}{!}{\includegraphics{./fig2.eps}}
\caption[]{Evolution of the \ha\ equivalent width. The circles represent
the data, and the solid lines are the best linear fit to the data points of the
corresponding interval. The red arrows mark the occurrence of an X-ray
outburst. The dash style simply denotes a weaker outburst.}
 \label{ewfit}
\end{figure}

\begin{table*}
\caption{Results of the linear fits to the different parts of the long-term
evolution of the \ew. See Fig.~\ref{ewfit} to identify the various regions.}
\label{linfit}
\center
\begin{tabular}{lccccc}
\hline
\hline
Region  &MJD range      &slope                  &Ordinate       &Corr. Coef.  &$N_{\rm obs}$  \\
\hline
\hline
Rise 1 (red)    &53544--55103   &$0.00598\pm0.00006$    &$-320\pm3$     &0.99   &30   \\
Rise 2 (blue)   &56220--58751   &$0.00410\pm0.00004$    &$-230\pm2$     &0.95   &67   \\
Decay 2 (green) &55103--56220   &$-0.00870\pm0.00011$   &$491\pm6$      &-0.92  &33   \\
Decay 3 (brown) &58751--59894   &$-0.00938\pm0.00010$   &$560\pm5$      &-0.98  &31   \\
\hline
\end{tabular}
\end{table*}
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig3.eps}
\caption[]{Representative profiles of the \ha\ line. }
\label{lineprof}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=10cm]{./fig4.eps}
\caption[]{Evolution of the \ha\ line parameters. {\em Top}: Separation between the blue and
red peaks (left axis) and the \ha\ EW (right axis).
{\em Bottom}: Evolution of the $V/R$ ratio.  }
\label{linevol}
\end{center}
\end{figure}

\subsection{\ha\ line profile variability}
\label{prof}

In contrast to what it is typically seen in other BeXBs, the shape of the \ha\ line
is rather symmetric.  While the strength of the line changes substantially over
time, as described in the previous section, the Be star in \gro\ always displays
a quite symmetric double-peak profile (Fig.~\ref{lineprof}). The ratio of the
peak intensities of the blue and red peaks is almost always $V/R\sim1$
(Fig.~\ref{linevol}, bottom panel), which it is rather unusual as most BeXBs
exhibit $V/R$ variability \citep{reig00,reig10b,yan16,alfonso17,zamanov20}.

The
separation between the blue and red peaks $\Delta V$ also evolves with time
 and is anticorrelated with \ew\ (Fig.~\ref{linevol}, top panel). At the optical
high states, the peak separation approaches the spectral resolution of our data
($\sim200$ km s$^{-1}$), and  the two peaks merge into a single-peak profile in
some spectra.

\subsection{Optical polarimetry}

The continuum polarization in Be stars is attributed to Thomson scattering of
unpolarized starlight in the disk
\citep{poeckert79,waters92,wood96,yudin01,halonen13a}.   In general, it is
difficult to compute the actual intrinsic polarization because most of the
observed polarization is expected to have an interstellar origin. For example,
using the relation between polarization and extinction $P_{\rm max,ISM} (\%) =
3.5 \times E(B-V)^{0.8}$ \citep{fosalba02}, we estimate that the maximum
contribution of the interstellar medium (ISM) to the measured optical
polarization toward \gro\ is $P_{\rm ISM}= 4.6$\%,  which is the observed
polarization, on average. However, the fact that  the polarization degree varies
by more than 1\% (from $\sim4$\% to $\sim5.2$\%) and that this variation
correlates with \ew\ (it increased as the \ew\ increased) supports the idea that
part of the observed polarization is intrinsic to the source and originates in
the disk. To estimate the contribution of the ISM, we observed three nearby
stars in the field of view of \gro. Table~\ref{polfsinf} lists the field stars
that we used, while Table~\ref{polfs} gives the polarimetric data for these
stars.

Figure~\ref{extdist} shows the change in extinction with distance and the
location of the field stars and the source. The data were taken from the  dust
maps by \citet{green19}. The figure indicates that at least three different
molecular clouds are in between the observer and the BeXB along the line of
sight. Clearly, as the distance increases, so does the polarization. The field
star fs2 is too close to Earth to represent a reliable measurement of the ISM
polarization at the source distance.  On the other hand, gfs2 displays  stronger
polarization than the source, which might indicate that the star is
intrinsically polarized and hence again only a poor representation of the ISM
polarization. We note that at very short and long distances, the computation of
the extinction is very uncertain due to the lack of sufficient number of stars.
If we use fs1 to estimate the polarization from the ISM (by vector subtraction,
the $q$ and $u$ Stokes parameters of the source and the field star), then we
obtain that the intrinsic polarization varied from 0.3\% to 1.3\%. The lowest
value of the polarization degree corresponds to the 2013 and 2022 observations,
when the \ew\ displayed minimum values.  The level of intrinsic polarization
that we measure agrees with the expected  polarization in an axisymmetric
circumstellar disk predicted by single-scattering plus attenuation models, which
is typically $\simless 2$\% \citep{waters92}. It also agrees with the measured polarization
in a sample of classical Be stars \citep{yudin01}.

\begin{table}
\caption{Selected field stars in the field of view of \gro\ used to derived the
ISM polarization. }
\label{polfsinf}
\begin{center}
\begin{tabular}{@~l@~c@~ccc}
\noalign{\smallskip}    \hline\noalign{\smallskip}
ID      &RA             &DEC           &Distance        &{\it Gaia} ID  \\
        &               &              &(pc)            &\\
\noalign{\smallskip}    \hline\noalign{\smallskip}
fs1     &314.72292      &41.77444       &$1894^{+67}_{-48}$     &2065652877361883648  \\
fs2     &314.69167      &41.76639       &$642^{+7}_{-6}$        &2065653598916381568  \\
gfs2    &314.77043      &41.76491       &$5787^{+1771}_{-1095}$ &2065652671203453696  \\
\noalign{\smallskip}    \hline\noalign{\smallskip}
\end{tabular}
\end{center}
\tablefoot{Distance from \citet{bailer-jones21}}
\end{table}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig5.eps}
\caption[]{Extinction as a function of distance (data from \citet{green19}).
The dots mark the assumed location of the stars on the extinction curve,
according to their distance. Distances are from \citet{bailer-jones21}. We also
indicate the polarization degree of each star (see Table~\ref{polfs}).}
\label{extdist}
\end{center}
\end{figure}

\subsection{Reddening and distance}

Owing to the presence of the circumstellar disk, the photometric magnitudes and
colors are affected by disk emission.  The overall effect is to make the star
appear redder, that is, to be seen under higher extinction than in the absence of
the disk. As a result, the estimate of the distance obtained from photometric
magnitudes that are thought to be contaminated by disk emission represents a lower
limit.  Some attempts have been made to quantify the contribution of the disk
to the photometric indices \citep{dachs88,fabregat90,riquelme12}. Here we take
advantage of the low optical states in \gro\ (Fig.~\ref{opt}). Because the
source of long-term optical variability is the disk, low optical states
correspond to instances in which the disk is almost absent. Therefore,  the low
optical states give us the opportunity to observe the underlying B star without
a major contribution from the disk.

To estimate the distance through the distance-modulus relation, $V-M_V-A_V=5
\log(d)-5$, the amount of interstellar extinction $A_{\rm V}=R \times E(B-V)$ to
the source has to be determined.  The most direct method for estimating $E(B-V)$ is
to use the calibrated color of the star according to the spectral type. The
color excess is defined as $E(B-V)= (B-V)_{\rm obs}-(B-V)_0$, where $(B-V)_{\rm
obs}$ is the observed color and $(B-V)_0$ is the intrinsic color of the star.
The expected color for an O9.5-B0V star is $(B-V)_0=-0.29\pm0.02$. This value is
the average of the calibrations from
\citet{johnson66}, \citet{fitzgerald70}, \citet{gutierrez79}, \citet{wegner94},
and \citet{pecaut13}. The error is
the standard deviation of the five values.   We applied this method to the three
low optical states independently and obtained $8.6\pm0.9$ kpc, $7.8\pm0.8$ kpc,
and $8.0\pm0.9$ kpc for the 2005--2006, 2012--2013, and 2022--2023 low states.
We note, however, that the average \ha\ EW was lowest during
the 2005--2006 period at \ew=--0.9 \AA, compared to \ew=--3.3 \AA\ and \ew=--1.3 \AA\
in the other two low states, respectively. Hence we take $8.5\pm0.8$ kpc as the
final value for the distance. This value compares to $8.7^{+1.1}_{-0.9}$ kpc
from {\it Gaia} \citep{bailer-jones21}. The average color excess is estimated to be
$E(B-V)=1.42\pm0.02$.

\section{Discussion}

BeXBs display variability on all timescales and at all wavelengths. The fastest
variability is detected in the X-rays (on the order of seconds) and corresponds
to the rotation of the neutron star, which manifests as pulsations. In the
optical band, the fastest timescales (on the order of hours or days) are
attributed to changes in the stellar photosphere (pulsation and rotation).
Long-term variations (on the order of months to years) are related to structural
changes in the disk. These long-term variations are the subject of the present
work.

\subsection{Disk formation and dissipation}

Be star disks are known to form and dissipate on timescales of years. The
disk emission increases with wavelength and becomes particularly significant in
the infrared band.
\gro\ went through two high (2008-2009 and 2018-2019) and three low optical
states (2004-2005, 2012-2013, and 2022-2023) during the course of our
observations. The transition from  low to high states and then back to low
states occurred smoothly over a period of several years. We attribute these
long-term changes to the formation and dissipation of the circumstellar disk.
As expected in the disks of Be stars, the amplitude of the variability in the
continuum increases as the effective wavelength increases: the difference in
magnitudes from maximum to minimum is $\Delta B=0.32$, $\Delta V=0.46$, $\Delta
R=0.56$, and $\Delta I=0.70$.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig6.eps}
\caption[]{Polarization angle and polarization degree of \gro\ for the period
2013-2022. The red arrow marks the occurrence of a type II outburst. Filled brown circles denote all the observations, and the open blue circles
represent the weighted mean calculated using the observations over one year. }
\label{pol}
\end{center}
\end{figure}

\subsubsection{Line and continuum emission}

In Sect.~\ref{res} we pointed out that the line emission variability timescales are faster than the continuum emission. The $BVRI$ magnitudes can remain
at a mean value with very small fluctuations during the low optical states. For
example, in the period MJD52800--54000, the mean and standard deviation of 14
measurements of the $V$ magnitude were 14.94 and 0.04 mag, respectively. For MJD
55800--56900 (six measurements), these values were 14.90 and 0.02 mag,
respectively. In contrast, \ew\ does not show a stable behavior like this and varies
more strongly. For the same periods, the average and standard deviation in
\ew\ was --1.8 and --0.8 \AA\ and --2.7 and --2.1 \AA, respectively.

Another interesting result is the fact that the peak in the \ew\ occurred about
400 days after the peak in the photometric magnitudes. Although the
observational gaps cannot be ignored, the excellent spectral coverage  clearly
indicates that the delay between the \ew\  and $R$-band maxima is real. These
delays imply that the continuum and line emission are anticorrelated for a certain
period of time (i.e., the time elapsed between the two peaks). This result has
been reported in many other BeXBs: Swift\,J0243.6+6124 \citep{liu22},
RX\,J0440.9+4431 \citep{yan16}, 1A\,0535+26 \citep{clark99,yan12b}, 4U\,0115+63
\citep{reig07b}, and SAX\,J2103.5+4545 \citep{camero14}. Interestingly, the
X-ray outbursts occur close to the \ew\ maxima.

This delay can be explained as due to the different locations in the disk in which
the line and continuum emission in Be stars originate. The outer parts of the
disk contribute most to the emission of the \ha\ line, whereas the origin of the
optical continuum resides in the inner parts. According to the viscous decretion
disk model \citep[see][for a review]{rivinius13a}, about 90--95\% of the
$V$-band flux comes from inside $1.8-2.5\,R_*$ \citep{haubois12}, while \ha\
emission reaches this percentage above $\simmore 6-10\, R_*$ \citep{carciofi11}.
Therefore, the fact that  the optical brightness has already started to decrease
by the time the \ew\ reached the maximum indicates that the dissipation of the
disk began from the inner parts of the disk. Disk dissipation occurs when the
mass ejection mechanism ceases. As the angular momentum supply stops, the inner
parts of the disk are reaccreted onto the star. However, viscosity still
transports angular momentum outward. Thus the outer parts of the disk continue
to expand, while the inner parts move inward \citep{haubois12,carciofi12}. Since
the flux in the $V$ (or $R$) band is mainly produced in the inner parts of the
disk, we expect the visual flux to start decreasing before the \ha\ line flux,
which mainly originates in the outer parts.

The long-term trend of the polarization degree also supports a disk
formation and dissipation cycle. Polarization in Be stars results from the
scattering of stellar radiation by electrons in the circumstellar disk.
The degree of the polarization provides information about the number of
scatterers. Therefore, we should expect a significant decrease in the
polarization degree during low optical states, as is indeed observed
(Fig~\ref{opt}). Since these low optical states represent instances when the
disk almost vanished, we can assume that the polarization degree roughly
corresponds with that of the ISM. With this assumption, we estimate that the ISM
polarization from the observer to \gro\ is $\sim4$\%.

Strictly speaking, none of low optical states corresponded to the complete loss
of the disk because the \ha\ line did not revert fully into absorption (see
Fig.~\ref{lineprof}).  A normal (i.e., non-\ha\ emitting) B0V star is expected to
have \ew$=+3.5$ \AA\  \citep{jaschek87}, while the lowest values of \ew\ during
the low optical states were around $-[0.3-1]$ \AA. Nevertheless, the low optical
states do represent highly debilitated disk phases.

\subsubsection{Connection with X-rays}

As we mentioned above, the \ew\ during the 2018-2019 maximum reached higher
values than during the 2008-2009 high state. Since \ew\ provides a measure of
the size of the disk, we conclude that the Be star developed a larger disk
during the 2018-2019 event. The accretion of a large amount of matter onto the
neutron star led to the type II outburst observed by the {\it NuSTAR} and Astrosat space
missions. In principle, it might be argued that since the disk was smaller during the
2008-2009 optical peak, the X-ray emission would be lower, as observed.   On the
other hand, given the similarities and regularity of the low and high optical
states (Figs.~\ref{opt} and \ref{ewfit}), the absence of a giant X-ray outburst
prior to the 2008--2009 optical maximum is surprising. The {\it Swift}/BAT hard X-ray
transient monitor detected an increase in 15--50 keV flux by a factor of two in
less than a week. A pointed observation with {\it Swift}/XRT was made on 9 May 2008.
The 0.3--10 keV flux was estimated to be $1.5 \times 10^{-10}$ erg cm $^{-2}$
s$^{-1}$, which at a distance of 9 kpc corresponds to a luminosity of $L_x=1.4
\times 10^{36}$ \ergs\ \citep{krimm08}. Although this luminosity is one order of
magnitude lower than the 2019 event, we cannot be sure that it corresponded to
the peak of the X-ray event. The 2008--2009 is reminiscent of a type I outburst,
and the lack of data can be attributed to the weakness of the detection due to
the large distance to the source.

One possible explanation for the different luminosity of the two events could be
a different geometry in a misaligned and tilted disk. Current models that
explain X-ray outbursts in BeXBs invoke warped, tilted, and misaligned disks.
Because disks in BeXBs are truncated \citep{reig97a,reig16}, the only way for
large amount of material to be transferred to the neutron star is in highly
asymmetric configurations in which a misaligned disk becomes warped and eccentric.
Highly distorted disks result in enhanced mass accretion when the neutron star
moves across the warped part \citep{martin11,martin14a,okazaki13}. Evidence of these warped disks has been reported for 1A\,0535+262 \citep{moritani13} based
on a spectral line analysis and for 4U\,0115 \citep{reig18b} based on variations in the polarization
angle. Although we observe a small change in the polarization angle
during the 2019 X-ray outburst (Fig.~\ref{pol}), it is not very significant.
Likewise, the symmetry of the \ha\ line profile argues against a highly
distorted disk. Therefore, it remains to be seen how a type I outburst can lead
to the almost complete destruction of the circumstellar disk.

We also noted in Sect~\ref{res} that the dissipation of the disk is faster than
its formation. The rate of change of the \ha\ equivalent width is about
1.8 \AA/yr during disk formation and about 3.3 \AA/yr during
disk dissipation (see Fig.~\ref{ewfit} and Table~\ref{linfit}). This result
contrasts with what it is observed in classical Be stars, in which the
timescales for disk growth are shorter than the timescales for disk dissipation
\citep{haubois12}. A crucial difference is the presence of the neutron star in
BeXBs. As the disk expands and reaches periastron distance, a substantial amount
of matter will be accreted onto the neutron star. Therefore, the dissipation of
the disk in BeXBs may occur faster owing to the interaction between the neutron
star and the disk.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig7.eps}
\caption[]{ ($(B-V)-V$ color magnitude diagram. Empty symbols correspond to
the disk formation phases, and filled symbols correspond to the dissipation
phases.}
\label{vbv}
\end{center}
\end{figure}

\subsection{Inclination angle}
\label{ang}

We showed in Sect.~\ref{prof} that the source displays a double-peak profile even
during the high optical state (see Fig~\ref{lineprof}). This result suggests that we
see the circumstellar disk at an intermediate inclination angle. Typically,
single-peak profiles are associated with high inclination angles (pole-on
stars) and shell profiles with high inclination angles (edge-on stars;
\citealt{slettebak79,hanuschik95,rivinius06,silaj14,sigut22}). An example of a low
inclination system is V\,0332+54, whose inclination angle is estimated to be
$i<20^{\circ}$ \citep{negueruela99,zhang05}, and a BeXB that exhibits shell
profiles is IGR\,J21343+4738 \citep{reig14a}.

The evolution of the photometric magnitudes and colors can also be used to
constrain the inclination angle. Figure~\ref{vbv} shows the color-magnitude
diagram of \gro. Empty symbols correspond to the disk formation phases, and
filled symbols correspond to the dissipation phases. The positive correlation
between the color index $(B-V)$ and the $V$ magnitude implies that the system is
viewed at an intermediate inclination angle, where the disk emission to the
optical colors and magnitudes contributes significantly. As the disk grows, the
system brightens ($V$ decreases), while the overall (disk plus star) emission
becomes redder ($(B-V)$ increases).  At very high inclination angles (equator-on
stars), the inner parts of the Be envelope partly block the stellar photosphere,
while the small projected area of the disk on the sky keeps the disk emission to
a minimum. Thus, stars viewed at very high inclination angles ($i > i_{\rm
crit}$) would show an inverse correlation \citep{harmanec83,harmanec00}.   The
value of the critical inclination angle is not known, but a rough estimate based
on available data suggests $i_{\rm crit} \sim 60-70^{\circ}$
\citep{hanuschik96a,sigut13}.
In summary, the correlation between the photometric colors and magnitudes and
the stable double peak profile of the \ha\ line in \gro\ suggest an inclination
angle in the range $30^{\circ}-60^{\circ}$.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig8.eps}
\caption[]{Peak separation as a function of \ew.  }
\label{deltav-ew}
\end{center}
\end{figure}

\subsection{Disk size and kinematics}
\label{velaw}

It was already proposed by \citet{struve31} that disks in Be stars are supported by rotation. The main evidence for this comes from the
correlation between the projected rotational velocity and the width and shape of
the emission line \citep{andrillat83,dachs86b,hanuschik88,hanuschik89,dachs92}.

In Sect~\ref{prof} we showed that the peak
separation $\Delta V$ and \ew\ are anticorrelated. This is the expected behavior in rotation-supported disks, in which the rotational velocity of the gas particles in the disk
varies as the inverse of the radius, $v_{\phi}\propto r^{-j}$  (e.g., for
Keplerian disks $j=1/2$).  In these disks, the peak separation roughly
corresponds to twice the projected rotational velocity of the particles in the
disk $\Delta V\approx 2 v_{\phi} \sin i$. As the radius of the disk increases
(i.e., as \ew\ increases), the rotational velocity decreases, and so does the peak
separation.  The index $j$ can be estimated from the relation between peak
separation and EW \citep{hanuschik88},

\begin{equation}
log(\Delta V/2 \sin i) = a \times \log(EW(H\alpha))+b
,\end{equation}

\noindent where the slope $a$ is related to the $j$ index as $a=j/2$. Figure \ref{deltav-ew} shows
this relation for \gro. When we restrict the fitting range to EWs between $1-8$ \AA, then we find $j=0.58\pm0.04$, which is not far from the
expected value of 0.5 for a Keplerian disk. Lower and higher values of the \ew\ were removed from the fit because below $\log(EW(H\alpha))\sim
0$ the disk is too small and has probably not reached a stable configuration,
while above 8 \AA, the $\Delta V$ saturates, which might be a spectral
resolution effect.

After showing that the particles in the disk of \gro\ follow a Keplerian
law, we can estimate the disk size. The disk velocity adopts the form
\citep[see. e.g.,][]{huang72}

\begin{equation}
\label{vphi}
v_\phi=v_0\left(\frac{r}{R_*}\right)^{-1/2}
,\end{equation}

\noindent where $R_*$ is the equatorial stellar radius, and $v_0$ is a initial
value of the disk velocity close to the stellar surface. In the limiting case,
$v_0$ would be the stellar rotational velocity. For simplicity, we
assumed $v_0=v_*$. By inverting eq.~(\ref{vphi}) and since
$v_{\phi} \approx (\Delta V/2 \sin i),$

\begin{equation}
\frac{R_d}{R_*} = \left(\frac{2 v_* \sin i}{\Delta V}\right)^2
,\end{equation}

\noindent where $r=R_d$ is the radius of the \ha\ emission region.
Figure~\ref{radew} shows the size of the circumstellar disk as a function of
time and \ew. As mentioned above, the leveling off at $\sim 50\,~{\rm R}_\odot$ marks the
minimum distance between the peaks that our spectrograph is capable of
discerning.

The Be star projected rotational velocity in \gro\ has been estimated by
\citet{kiziloglu07b} to be in the range $240-310$ km s$^{-1}$. Here we assumed $v_*
\sin i = 275$ km s$^{-1}$ , which corresponds to half the maximum peak separation.
We note, however, that \ha\ double-peak separation exceeding $2 v_* \sin i$ has
been measured in Be stars during phases of very weak \ha\ emission
\citep{dachs92}.

The largest disk radius is estimated to be $R_{\rm disk,max}\sim5-6 \, R_*$ or
$60\, ~{\rm R}_\odot$, assuming a stellar radius of an O9.5-B0e IV-V star of $10\, ~{\rm R}_\odot$
\citep{martins05}.  From Kepler's third law, we can estimate the orbital
separation $a\sim 160 \, ~{\rm R}_\odot$ for $P_{\rm orb}=55$ days and typical values of
the mass $M_*=18\, ~{\rm M}_\odot$. If we now assume that in order to produce an X-ray
outburst, the neutron star must interact with the disk, that is, the disk radius
must be similar to the periastron distance ($a(1-e)$), then the eccentricity of
the system should be $e\sim 0.6$.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig9.eps}
\caption[]{Disk radius as a function of \ew. }
\label{radew}
\end{center}
\end{figure}

\section{Conclusion}

We have performed the most detailed study of the long-term optical variability
of the BeXB \gro.  The system displays correlated variability on timescales of
years in all the optical observables. The source goes through optical high and
low states. The high optical state is characterized by large \ew\, bright
continuum emission, and significant polarization. We identify this state with a
well-developed decretion disk. The smooth increases and decreases of the optical
emission are interpreted as phases of formation and dissipation of the
circumstellar disk. The optical maxima coincide with the occurrence of X-ray
outbursts, after which all optical indicators decreased, leading to the
dissipation of the disk.  However, a disk-loss episode never occurred.  The
almost permanent double-peak line profile indicates that the systems is seen at
an intermediate inclination angle, whereas the absence of $V/R$ variability is
attributed to a relatively small and stable disk. We did not find evidence for a
warped disk during the 2019 type II X-ray outburst.

\begin{acknowledgements}

Skinakas Observatory is run by the University of Crete and the Foundation for
Research and Technology-Hellas.   The Liverpool Telescope is owned and operated
by the Astrophysics Research Institute of Liverpool John Moores University. The
Starlink software is currently supported by the East Asian Observatory. IRAF is
distributed by the National Optical Astronomy Observatories, which is operated
by the Association of Universities for Research in Astronomy, Inc. (AURA) under
cooperative agreement with the National Science Foundation. This research has
made use of the SIMBAD database, operated at CDS, Strasbourg, France and of
NASA's Astrophysics Data System operated by the Smithsonian Astrophysical
Observatory.

\end{acknowledgements}

\bibliographystyle{aa}


\section*{\refname}}

\usepackage[usenames]{color}
\newcommand{\B}[1] {{\color{blue} #1}}
\newcommand{\R}[1] {{\color{red} #1}}
\newcommand{\G}[1] {{\color{green} #1}}

\begin{document}

	\title{Supporting Information for: \\
 Transient 2D IR spectroscopy and multiscale simulations reveal vibrational couplings in the Cyanobacteriochrome Slr1393-g3}

\author{David Buhrke$^{1,2,*}$, Yigal Lahav$^{3,4}$, Aditya Rao$^{3}$, Jeannette Ruf$^{1}$, Igor Schapiro$^{3}$ and Peter Hamm$^{1}$\\
\textit{$^1$Department of Chemistry, University of Zurich, Zurich, Switzerland\\
$^2$current adress: Institute of Biology, Humboldt University Berlin, Germany\\
$^3$Fritz Haber Center for Molecular Dynamics Research Institute of Chemistry, The Hebrew University of Jerusalem, Israel\\
$^4$ MIGAL - Galilee Research Institute, Kiryat Shmona, Israel\\}
$^*$corresponding author: david.buhrke@hu-berlin.de
}

\maketitle

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{S1.png}
	\caption{Displacement vectors of the calculated vibrational modes in the  Pg state (left, green) and the Pr state (right, pink).}
 \label{intro}
\end{figure}

\begin{table}
\caption{Calculated anharmonicity of the overtones (diagonal, defined as the difference between overtone and fundamental transition frequency, $\omega_{1,2}-\omega_{0,1}$) and combination bands (off-diagonal, defined as the difference between combination mode and fundamental transition frequency, $\omega_{01,11}-\omega_{00,10}$), in cm$^{-1}$ of the relevant vibrational modes in the Pg state.}

\begin{tabular}{ |p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|  }
\hline
  & CC$_\mathrm{D}$  & CC$_\mathrm{BC}$ & CC$_\mathrm{CD}$ & CC$_\mathrm{AB}$ & CO$_\mathrm{D}$  &  CO$_\mathrm{A}$ \\
 \hline
CC$_\mathrm{D}$  & -12.04 & -1.28 & -1.94 & -0.27 & -1.67 & 0.02 \\
\hline
CC$_\mathrm{BC}$ &  & -0.85 & -0.17 & -4.78 & -0.35 & -0.64 \\
\hline
CC$_\mathrm{CD}$ &  &  & -11.01 & -0.12 & -1.95 & 0.00 \\
\hline
CC$_\mathrm{AB}$ &  &  &  & -4.50 & -0.19 & -1.00 \\
\hline
CO$_\mathrm{D}$  &  &  &  &  & -17.13 & 0.03 \\
\hline
CO$_\mathrm{A}$  &  &  &  &  &  & -17.23 \\
  \hline
\end{tabular}
\label{TablePg}
\end{table}

\begin{table}
\caption{Calculated anharmonicity of the overtones (diagonal, defined as the difference between overtone and fundamental transition frequency, $\omega_{1,2}-\omega_{0,1}$) and combination bands (off-diagonal, defined as the difference between combination mode and fundamental transition frequency, $\omega_{01,11}-\omega_{00,10}$), in cm$^{-1}$ of the relevant vibrational modes in the Pr state.}

\begin{tabular}{ |p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|  }
\hline
  & CC$_\mathrm{BC}$  & CC$_\mathrm{D}$ & CC$_\mathrm{AB}$ &CC$_\mathrm{CD}$ & CO$_\mathrm{D}$  &  CO$_\mathrm{A}$ \\
 \hline
CC$_\mathrm{BC}$  & 0.46 & -2.14 & -0.56 & -0.63 & -0.57  & -0.37 \\
\hline
CC$_\mathrm{D}$  &  & -7.56 & -1.25 & -3.95 & -1.42 & 0.00 \\
\hline
CC$_\mathrm{AB}$ &  &  & -4.98 & -2.46 & -0.48 & -2.16 \\
\hline
CC$_\mathrm{CD}$ &  &  &  & -6.08 & -0.66 & -0.36 \\
\hline
CO$_\mathrm{D}$  &  &  &  &  & -16.84 & 0.06 \\
\hline
CO$_\mathrm{A}$  &  &  &  &  &  & -15.67 \\
  \hline
\end{tabular}
\label{TablePr}
\end{table}

\date{\today}

\end{document}

\section*{\refname}}

\usepackage[usenames]{color}
\newcommand{\B}[1] {{\color{blue} #1}}
\newcommand{\R}[1] {{\color{red} #1}}
\newcommand{\G}[1] {{\color{green} #1}}

\begin{document}

	\title{Transient 2D IR spectroscopy and multiscale simulations reveal vibrational couplings in the Cyanobacteriochrome Slr1393-g3}

\author{David Buhrke$^{1,2,*}$, Yigal Lahav$^{3,4}$, Aditya Rao$^{3}$, Jeannette Ruf$^{1}$, Igor Schapiro$^{3}$ and Peter Hamm$^{1}$\\
\textit{$^1$Department of Chemistry, University of Zurich, Zurich, Switzerland\\
$^2$current adress: Institute of Biology, Humboldt University Berlin, Germany\\
$^3$Fritz Haber Center for Molecular Dynamics Research Institute of Chemistry, The Hebrew University of Jerusalem, Israel\\
$^4$ MIGAL - Galilee Research Institute, Kiryat Shmona, Israel\\}
$^*$corresponding author: david.buhrke@hu-berlin.de
}

\date{\today}

\begin{abstract}
Cyanobacteriochromes are bi-stable photoreceptor proteins with desirable photochemical properties for biotechnological applications such as optogenetics or fluorescence microscopy. Here, we investigated Slr1393-g3, a cyanobacteriochrome that reversibly photo-switches between a red-absorbing (Pr) and green-absorbing (Pg) form. We applied advanced IR spectroscopic methods to track the sequence of intermediates during the photocycle over many orders in magnitude in time. In the conversion from Pg to Pr, we have revealed a new intermediate which precedes the Pr formation by using transient IR spectroscopy. In addition, stationary and transient 2D~IR experiments measured the vibrational couplings between different groups of the chromophore and the protein during these intermediate states. Anharmonic QM/MM calculations predict spectra in close-to-quantitative agreement with experimental 2D~IR spectra of the initial and the final state of the photocycle. They facilitate the assignment of the IR spectra and provide an atomistic insight into the coupling mechanism. This serves as a basis for the interpretation of the spectroscopic results and suggests structural changes of the intermediates along the photocycle.
\end{abstract}

\maketitle

\section{Introduction}

\noindent Cyanobacteriochromes (CBCRs) are bi-stable cyanobacterial photoreceptors\cite{Rockwell2010, Fushimi2019} with emerging applications in biotechnology, such as in super resolution microscopy\cite{Oliinyk2019} or optogenetics.\cite{Blain-Hartung2018} CBCRs typically consist of multiple photosensory modules that each incorporate an open-chain tetrapyrrole chromophore such as phycocyanobilin (PCB). In the present study, we investigated Slr1393-g3 (from here on denoted Slr-g3 for simplicity), a photosensory domain from a red/green CBCR that incorporates PCB and converts reversibly between a red-absorbing (Pr) and a green-absorbing (Pg) parent state.\cite{Chen2012, Slavov2015, Xu2014, Xu2020}
Both states are modulated by geometric changes of the PCB chromophore that are stabilized by the protein environment.
Figure \ref{intro} shows the structure of PCB inside the binding pocket of Slr-g3, which adopts a highly distorted geometry in Pg, thereby shortening the effective conjugation length compared to Pr.\cite{Xu2020, Wiebeler2018} The exact course of events that leads to these geometric changes is not fully established, hence a detailed mechanistic understanding of the photocycle reactions in red/green CBCRs is still lacking. Time-resolved spectroscopy studies indicated that these proteins show large shifts in their visible absorption properties on the \textmu\text{s}-ms time scale after photoexcitation\cite{Xu2014, Jenkins2019, Kirpich2018}, but it is up to now largely unclear how exactly PCB changes its structure as a function of time and how this process is driven by the protein environment.

\begin{figure*}[htbp]
	\centering
	\includegraphics[width=0.75\linewidth]{QM_intro.png}
	\caption{The QM/MM partitioning and optimized structures of Slr-g3 in the Pg state (left panel) and Pr state (right panel). The chromophore is comprised of 86 atoms (shown in licorice representation). The pyrrole hydrogens were deuterated (shown in dark gray) and positions of the propionate atoms were fixed (shown in light gray) during the calculation of vibrational frequencies. The nomenclature of the four pyrrole rings \textit{A-D} and localisation of the relevant normal modes are indicated in the right panel.
 }
 \label{intro}
\end{figure*}

Infrared (IR) spectroscopy is a powerful technique to study exactly such protein reaction dynamics, i.e. the light-induced cycles of photoreceptors or catalytic cycles of enzymes.\cite{Kottke2017,Lorenz-fonfria} One of the most popular observables in the mid-IR is the so-called amide I region between 1600 and 1700~cm$^{-1}$, that is dominated by signals from the peptide bond C=O stretching vibrations (called amide I modes).\cite{Barth2002} Additionally, co-factors such as PCB reveal structure-sensitive vibrational modes in the same spectral region, resulting in small IR signals superimposed on the overwhelming amide I band. In photoreceptors, these are e.g. the C=C, C=O or C=N stretching vibrations of retinal, tetrapyrrole or other co-factors that can be identified in light-dependent difference-IR experiments. The same vibrational modes can be measured in resonance Raman (RR) experiments, where they are resonance enhanced over the amide I background. Thus,  RR and IR can be used for cross-validation and provide assignments for difference-IR experiments.\cite{Buhrke2019a, Lorenz-fonfria} The amide~I region of different proteins was investigated in detail not only by linear (1D) but also two dimensional (2D-)IR spectroscopy, which offers additional information such as spectral diffusion, vibrational couplings or anharmonicities.\cite{GANIM2007} Generally, the analysis of the amide I region is tedious due to a strong coupling of the broad amide I modes, and due to spectral congestion, which makes it extremely difficult to extract local information even for small proteins and peptides.\cite{Hamm1998, Strasfeld2009, GANIM2007} Due to this obstacle, most applications of 2D~IR spectroscopy on protein samples require carefully designed experimental conditions such as isotope labeling and very small systems,\cite{Bredenbeck2003,Shim2009} or the introduction of IR labels that absorb outside the amide vibrations, typically in the so-called ``transparent window'' between 1800-2300~cm$^{-1}$.\cite{Adhikary2017,Thielges2021} Although 2D~IR spectroscopy has been thoroughly developed over the last two decades, the technique is still non-standard and requires highly specialized labs and on top of that very expensive equipment such as femtosecond laser systems and optical parametric amplifiers. Due to a combination of these challenges, 2D~IR spectroscopy has so far been applied  only sparsely to investigate the light-induced reactions of photoreceptor proteins, e.g. on Bacteriorhodopsin\cite{Andresen2009, Hamm2021}, Photoactive Yellow Protein\cite{Schmidt-Engler2020}, a CBCR\cite{Ruf2021}, and a Phytochrome.\cite{Buhrke2022}

Here, we tracked the IR response of Slr-g3 after excitation of Pg in the \textmu\text{s}-ms time range with transient IR spectroscopy; we will call this rather conventional method transient 1D~IR spectroscopy from now on. Further, we show that the transformation of the chromophore geometry is accompanied by changes in the coupling between certain vibrational normal modes that show up as difference cross-peaks in a stationary difference-2D~IR experiment. We then combined these two approaches and investigated the photorecation with transient 2D~IR spectroscopy in the \textmu\text{s}-ms time range.\cite{Hamm2021} Through the combination of these three related IR methods with quantum mechanics/molecular mechanics (QM/MM) simulations and anharmonic frequency calculations, we were able to develop a solid mechanistic model for the photocycle of Slr-g3.

\section{Materials and Methods}

\subsection{Protein expression and purification}

Slr-g3 was expressed in \textit{E. coli} with an \textit{in vivo} chromophore assembly protocol and purified as described previously.\cite{Buhrke2020b} For all IR experiments, the samples were dissolved in 50~mM Tris, 200~mM NaCl D$_2$O buffer (pD=7.8) and had a concentration of ca. 2~mM as determined by their UV absorption at 280~nm.

\subsection{IR Spectroscopy}

The setup used for all IR experiments was described elsewhere in detail.\cite{Hamm2021} Briefly, a commercial 100~kHz laser system (Amplitude) and optical parameteric amplifier (FastLite) were used to generate fs mid-IR pulses. The IR pulses were split into pump, probe and reference beams, the latter two of which were passed through the sample cuvette separated by roughly 1~mm, and then imaged onto the entrance slit of a spectrograph that was equipped with a 2x32-element MCT array detection and single-shot ADC electronics at 100~kHz.\cite{Farrell2020} In most experiments, the enhanced referencing scheme introduced in Ref.~\onlinecite{Feng2017} has been applied.

For transient 1D~IR spectroscopy, the IR pump pulses were blocked, while green actinic pulses were derived from a Nd:YAG laser (AO-S-532, CNI, Changchun, China, pulse duration ca. 40~ns) that was electronically synchronized to the 100 kHz probe light. Transient 1D~IR spectra were recorded in a vis-pump-IR-multiprobe scheme.\cite{Greetham2016,Hamm2021,Jankovic2021b,Buhrke2021} The green pump pulse energy was 17~\textmu\text{J} and the polarization was set to magic angle relative the probe beam to suppress orientational effects.

For stationary 2D~IR spectroscopy, the actinic pump pulse was blocked, while a pair of IR pump pulses generated in a pulse shaper (PhaseTech) were imaged onto the sample.\cite{Middleton2010} We used 4-state phase cycling to suppress scattering and strong undersampling in a rotating frame to minimize the number of laser shots (down to 168) needed to measure a complete 2D~IR spectrum. The waiting time between pump pulses and probe pulse was set to 200~fs.

Finally, for transient 2D~IR spectroscopy, both the IR-pump pulse pair as well as the actinic pump were imaged onto the sample. We used an interleaved sampling of 2D~IR and actinic pump delays, as introduced in Ref.~\onlinecite{Hamm2021}, which allowed us to generate a sequence of transient 2D~IR  spectra by proper reshuffling the data, each separated by 10~\textmu\text{s}. The transient 2D~IR  spectra were binned on a logarithmic time-scale with 10 time points per decade, rendering the signal-noise ratio at later delay times better than at earlier times since more transient 2D~IR  spectra enter a bin.

For the transient 1D~IR experiment, the repetition rate of the actinic pump laser, which determines the measurable time window, has been 15~Hz. It has been increased to 30~Hz for the transient 2D~IR  measurements, just enough to cover the complete reaction cycle, again in order to minimize the number of laser shots needed to measure a complete sequence of 2D~IR spectra (ca. 600000).

The sample was contained in a closed-cycle flow system consisting of a peristaltic pump and a sample flow cell with solenoid-driven micro valves developed for precise sample exchange synchronized to the laser system.\cite{Buhrke2021} The flow cycle contained an additional cell between the peristaltic pump and the pressure reservoir that allowed to prepare samples in the Pg state by illumination with a red laser diode (HL6750MG, Thorlabs) before entering the measurement cell (path length 50~\textmu\text{m}).

\subsection{QM/MM calculations}

The model preparation based on the crystal structure of Slr1393g3\cite{5dfx} can be found in detail elsewhere.\cite{Wiebeler2019,Wiebeler2019_2,Rao2021} Briefly, the structures of both Pr and Pg forms (PDB IDs: 5DFX and 5M82, respectively) were protonated and embedded in a TIP3P\cite{tip3p} water box, followed by MM optimization keeping the chromophore restrained to its crystal structure geometry. In this work, starting with the MM optimized structures, all solvent atoms distant more than 3.5~$\mathrm{\AA}$ of the chromophore were removed from the structure and a preliminary QM/MM optimization was performed at XTB2\cite{Grimme19}:ff14SB\cite{Simmerling15} level of theory using the ORCA\cite{ORCA} software package. The QM region contains 86 atoms (Fig. \ref{intro}), which includes the chromophore and the side chain of its binding cysteine residue. A hydrogen link atom is used to saturate the truncated bond between the C$_{\alpha}$ and C$_{\beta}$ carbons of the cysteine sidechain. A second QM/MM optimization was carried out at the B3LYP\cite{Becke92}/6-31G*\cite{Pople98}:ff14SB level of theory with Grimme's D3 dispersion correction with Becke-Johnson damping\cite{Grimme11} using the Gaussian16\cite{Gaussian16} software package. During this step, all MM atoms were kept fixed to their spatial positions, and the PCB pyrrole hydrogens were replaced by deuterium atoms. The optimized structures were subject to IR spectra calculation using the harmonic approximation. In this step the propionate groups of the chromophore were fixed. For calculating the 2D IR spectrum we followed the procedure by Mukamel et al.\cite{Mukamel2003,Mukamel2003_2}, focusing on the harmonic modes between 1600 and 1800~cm$^{-1}$. These include the two C=O stretches of rings A and D, the three C=C stretches of the methine bridges, and the C=C stretch inside ring D (right panel of Fig. \ref{intro}). These modes were subject to anharmonic frequencies, overtones and combination bands, calculation using second order perturbation theory\cite{Barone2004,Barone2005}, as implemented in the Gaussian16 package. It was validated by Geva et. al.\cite{Geva09} against full diagonalization of the vibronic Hamiltonian and was within good agreement. The final 2D IR signal was calculated by non-linear response theory\cite{hamm11} using the anharmonic frequencies and the harmonic IR intensities. The homogeneous dephasing time T$_2$ was taken as 1.0~ps.

\begin{figure*}[htbp]
	\centering
	\includegraphics[width=0.9\linewidth]{2DIR.png}
	\caption{Stationary 1D~IR and 2D~IR spectroscopy of Slr-g3 in the amide I region. (A,B) 1D~IR spectra in the Pg and Pr states, respectively, and (C) the ``Pr-minus-Pg'' difference 1D IR spectrum. The blue line in panel (C) compares the staionary difference spectrum with the late time transient difference spectrum measured 11~ms after excitation. Panels (D)-(F) show the corresponding 2D~IR spectra, where blue colors indicate the negative signal for the 0-1 transition of the corresponding oscillators, and red colors positive signals from the anharmonically shifted 1-2-transitions. Calculated 1D~IR spectra of PCB in the Pg and Pr states, as well as the ``Pr-minus-Pg'' difference spectrum are shown in panels (G)-(I) and the corresponding 2D~IR spectra in panels (J)-(L).}
	\label{2DIR}
\end{figure*}

\section{Results}

\subsection{Stationary difference spectroscopy}

\begin{table*}
\caption{Experimental and calculated anharmonic frequencies and intensities of relevant vibrational modes in the Pg/Pr states.}
\begin{tabular}{ |p{2.5cm}||p{3cm}|p{3.5cm}|p{4cm}|p{6cm}|  }
\hline
Mode & Exp. freq. (cm$^{-1}$)  & Calc. freq. (cm$^{-1}$) & Calc. IR int. (km mol$^{-1}$) \\
 \hline
CO$_\mathrm{A}$   & n.d./1720  &  1762/1724 & 356/310\\
CO$_\mathrm{D}$   & 1700/1675  &  1697/1670 & 714/1117\\
CC$_\mathrm{AB}$  & 1630/1627  &  1635/1619 & 647/846\\
CC$_\mathrm{CD}$  & n.d./n.d.  &  1617/1632 & 203/192\\
CC$_\mathrm{BC}$  & 1603/1600  &  1612/1595 & 2301/1637\\
CC$_\mathrm{D}$  & n.d./n.d.   &  1605/1594 & 12/24 \\
 \hline
\end{tabular}
\label{Table}
\end{table*}

The 1D IR absorption spectrum of Slr-g3 is dominated by amide I vibrations in the form of an asymmetric, broad and unstructured band in the region between 1600 and 1700~cm$^{-1}$. The spectra in the Pg and Pr states appear to be identical because the spectral changes associated with the photoconversion are too small to see them on the scale of Fig.~\ref{2DIR}A and B. The small changes associated with the photoconversion can be extracted by calculating a  ``Pr-minus-Pg'' difference 1D~IR spectrum, see Fig.~\ref{2DIR}C. Here, negative features correspond to loss of absorption due to Pg bleach, while positive bands indicate an absorption gain in Pr.

The 2D~IR spectra of Slr-g3 in the Pr and the Pg states are  indistinguishable as well in the representation of Fig~\ref{2DIR}D and E. All peaks in a 2D~IR spectrum appear as pairs with a negative (blue) lobe at higher probe frequencies and a positive (red) lobe at lower probe frequencies. Due to the quadratic dependence of the 2D~IR signal on the absorption cross section and the resulting increase in spectral resolving power, two peaks at 1628~cm$^{-1}$ and 1645~cm$^{-1}$ can be distinguished on the diagonal that only appear as slight asymmetry in the 1D~IR absorption spectra. These main amide I peaks are tentatively assigned to contributions of $\alpha$-helical and $\beta$-sheet structures of the protein, respectively.\cite{Barth2007}

Similar to the difference 1D~IR spectrum, a ``Pr-minus-Pg'' difference 2D~IR spectrum was calculated, which bears information about light-induced structural changes (Fig. \ref{2DIR}F). All features in the difference 1D~IR spectrum have counterparts on the diagonal of the difference 2D~IR spectrum that correspond to the change in IR absorbance at the respective frequencies caused by the photoconversion of the protein. However, each diagonal peak consist of two contributions, one on the diagonal corresponding to the 0-1-transition of the corresponding vibrational mode and a second one below the diagonal with opposite sign due to the anharmonically shifted 1-2-transition (excited state absorption).\cite{hamm11} Depending on the signs of the various bands in the difference 1D~IR spectrum, the signs of the diagonal peaks vary.

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.75\linewidth]{Raman.png}
	\caption{Vertical cuts through the experimental difference 2D~IR spectrum at 1600~cm$^{-1}$ (magenta) and 1627~cm$^{-1}$ (orange). The respective dotted lines are taken from a separate ``Pg-minus-Pr'' control experiment for reversibility. Red line: RR spectrum of Slr-g3 in the Pr state (reproduced from ref.~\onlinecite{Buhrke2020}).}
	\label{Raman}
\end{figure}

In addition, cross peaks are observed in the difference 2D~IR spectrum, labelled  \textbf{III},  \textbf{IV} and \textbf{V} in Fig. \ref{2DIR}F. The signs of all these cross peaks indicate they originate from Pr. In order to verify these cross peaks despite their small size (they are more than a factor 100 smaller than the 2D-IR spectra of Pg and Pr), the sample was converted back and forth between both states multiple times, and the inverted ``Pg-minus-Pr'' difference spectrum was also calculated. The reversibility and relative size of the cross-peaks can best be judged in vertical cuts at $\omega_{pump}$ at their maxima shown in Fig.~\ref{Raman}.  Here, the solid magenta and orange lines represent cuts through the ``Pr-minus-Pg'' spectrum while dashed lines correspond to an independently measured, reversed ``Pr-minus-Pg'' experiment.

We also compare these cuts in Fig.~\ref{Raman} with the experimental resonance Raman (RR) spectrum of Slr-g3 in the Pr state that we reproduced from Ref.~\onlinecite{Buhrke2020}. The frequencies of the cross- and diagonal peaks completely coincident with the main peaks of the RR spectrum, confirming that these signals originate from the conjugated C=C system. In Fig.~\ref{2DIR}F, the difference cross-peaks are clearly visible only on one side of the main diagonal with  $\omega_{probe}>\omega_{pump}$, while they are masked by the stronger excited state absorption contributions of the diagonal signals on the other side.  Nevertheless, in Fig.~\ref{Raman}, a superimposed cross-peak can be identified as a kink at 1600~cm$^{-1}$ in the vertical cut at 1627~cm$^{-1}$ (orange line, labeled \textbf{III}).

We anticipate that most features in the difference spectra originate from the pronounced geometric changes of the PCB chromophore, while some may also originate from the changes of protein structure that have impact on the amide I spectrum.\cite{Buhrke2020b} These signals generally overlap with each other. To distinguish these contributions, and to facilitate an assignment, we performed anharmonic QM/MM calculations where the protein environment was fixed and only vibrational modes from the chromophore were considered. Six vibrational modes of PCB were found in the investigated frequency window between 1580 and 1750~cm$^{-1}$, the two, mostly localized C=O stretching modes of the ring A and D carbonyls as well as C=C stretching modes of the conjugated system. The latter are  dominated by stretching motions of the three methine bridges (A-B, B-C and C-D), and are denoted accordingly as CC$_\mathrm{XY}$ in the following. One mode inside ring D (CC$_\mathrm{D}$, see Table~\ref{Table}) was identified as well, but its intensity is very weak, hence this mode will not be discussed here. The calculated atomic displacement vectors of all identified modes in Pg and Pr are given in Fig.S1.

The calculated frequencies and  IR intensities are listed in Table~\ref{Table} for both states Pg/Pr, the calculated anharmonicities are listed in Tables S1 and S2. Fig.~\ref{2DIR}G and H plot the resulting 1D~IR spectra, which are dominated by 3 peaks (i.e., CC$_\mathrm{BC}$ and CC$_\mathrm{AB}$ as well as CO$_\mathrm{D}$). CC$_\mathrm{CD}$ is very weak and visible only in Pr as a shoulder, while CO$_\mathrm{A}$ is outside the spectral window of Fig.~\ref{2DIR}. With the exception of  CC$_\mathrm{CD}$, all modes reveal a frequency downshift upon the Pg-to-Pr transition.

Based on the QM/MM results, and in combination of the Raman spectrum of Pr shown in Fig.~\ref{Raman},\cite{Buhrke2020} we can assign the experimental peaks of the Pr state, see Table~\ref{Table} (we will return to the Pg state later). This assignment is based on the energy ordering and relative intensities of the  various modes. It is interesting to note that among the bridging modes,  CC$_\mathrm{BC}$ has the strongest IR intensities while CC$_\mathrm{AB}$ has the stronger Raman intensity (Fig.~\ref{Raman}). In accordance with the QM/MM results, this implies that both modes are delocalized among each other to a certain degree, with one being more of a symmetric character and the other of a asymmetric character. Based on the results of the calculations, we conclude that CC$_\mathrm{CD}$ is too weak to be identified experimentally.

Each of the calculated resonances in the 1D~IR spectra shows up with the characteristic doublet on the the diagonal of the corresponding 2D~IR spectrum of Pg (Fig.~\ref{2DIR}J). In addition, strong cross peaks \textbf{I} (between CO$_\mathrm{D}$ and the CC modes) and \textbf{II} (between CC$_\mathrm{BC}$ and CC$_\mathrm{AB}$) reveal that the chromophore modes are coupled among each other.

The 2D~IR spectrum of Pr is better resolved (Fig.~\ref{2DIR}K), since the C=C resonances are spread over a wider frequency range. In particular, a cross peak \textbf{III}  can be identified, connecting the strongly coupled CC$_\mathrm{BC}$ and CC$_\mathrm{AB}$ modes. Furthermore,  CO$_\mathrm{D}$ now produces two distinct cross peaks with the CC$_\mathrm{BC}$ (\textbf{IV}) and CC$_\mathrm{AB}$ (\textbf{V}). These spectral changes give rise to a specific fingerprint pattern in the calculated difference spectrum (Fig. \ref{2DIR}L), reproducing the main features of the experimental difference spectrum surprisingly well (Fig. \ref{2DIR}F). This is in particular true for the cross peaks \textbf{III}, \textbf{IV} and \textbf{V} that connect the major signals from the C=C region with  CO$_\mathrm{D}$ in the Pr state. The  cross peak \textbf{I} with opposite sign, origination Pg, is not visible in the experimental spectrum, either because it is too weak or outside the measured frequency window.

One curiosity is the diagonal peak of CC$_\mathrm{BC}$, whose sign is inverted which is indicative for a positive anharmonicity with the 1-2 frequency being larger than the 0-1 frequency (Fig.~\ref{2DIR}K as well as Table SII in SI). A positive anharmonicity appears when the corresponding potential is steeper than harmonic (more particle-in-a-box like). With the exception of a proton vibrating between two heavier partners,\cite{Dahms2017, Fournier2018, Dereka2021} such as in the Zundel ion, a positive anharmonicity has not been observed experimentally. Here, it is the vibration of the bridging carbon between two heavier pyrrole rings that is predicted to have a positive anharmonicity, but only in the more rigid Pr state, and in the middle of the molecule (i.e., only for CC$_\mathrm{BC}$). The signature of a  positive anharmonicity is also visible in the calculated Pr-minus-Pg difference spectrum   (Fig.~\ref{2DIR}L), but it is not obvious in the experimental counterpart (Fig.~\ref{2DIR}F).

\begin{figure}[t]
	\centering
	\includegraphics[width=\linewidth]{T1DIR_V3.png}
	\caption{Transient 1D~IR spectroscopy of the Pg $\rightarrow$ Pr reaction from 10~\textmu\text{s} to 30~ms. (A) Data represented as a contour plot; blue colors indicate negative bands (bleaches) and red colors positive bands. (B) Dynamical content. (C) Three representative kinetic traces with lifetime analysis fits (black lines).}
	\label{T1DIR}
\end{figure}

\subsection{Transient 1D~IR and 2D~IR spectroscopy}

The stationary 1D~IR and 2D~IR difference spectra of Figs.~\ref{2DIR}C and F report on changes in the geometric and electronic structure between the two stable states of the protein. The time-resolved counterparts of these spectroscopies can provide the same information for transient intermediates during the course of the photocycle. In a recent publication,\cite{Buhrke2020b} we presented transient 1D~IR data of Slr-g3 in the ps to \textmu\text{s} time range obtained with a pump-probe setup based on two electronically synchronized femtosecond laser systems.\cite{Bredenbeck2004} Here, we add new data in the \textmu\text{s}-ms regime acquired with a conceptually different vis-pump-IR-multiprobe spectrometer, operating at repetition rates of 15~Hz (vis pump) and 100~kHz (IR probe, see Materials and Methods).\cite{Hamm2021,Jankovic2021b,Buhrke2021}
The new data reveal two further reaction steps that can be identified directly by visual inspection of the contour plot of Fig.~\ref{T1DIR}A, or in the three selected traces shown in Fig.~\ref{T1DIR}C.  Time constants for these processes were obtained by lifetime analysis with a maximum entropy method,\cite{Kumar2001,Lorenz-Fonfria:06} similar to the analysis that was applied for the first part of the data in our recent paper.\cite{Buhrke2020b} The resulting dynamical content $D$ is plotted in Fig.~\ref{T1DIR}B, which reveals two clear maxima at 250~\textmu\text{s} and 2.5~ms associated with the time constants of these processes.
The properly scaled late time (i.e., 11~ms) transient 1D~IR spectrum is overlaid in Fig.~\ref{2DIR}C (blue line) with the stationary 1D~IR difference spectrum (black line). The perfect agreement of both spectra indicates that the photoreaction is completed after that time and Pg is formed.

 The micro-to-milisecond kinetics of Slr-g3 were up to now only studied by transient visible spectroscopy at single wavelengths, which yielded only one time constant of 1.1~ms in this time window.\cite{Xu2014} This indicates that the two processes we see in the IR could not be distinguished in the visible experiment and a mix of the two time constants was observed. The same would be true if only single traces in the IR, e.g. at 1650 or 1725~cm$^{-1}$ were analyzed (see Fig.~\ref{T1DIR}A), where the two processes have the same sign.
 However, other traces, e.g. at 1610~cm$^{-1}$, show two clearly distinguishable processes with opposite sign and an absorption maximum around 1~ms (marked with an asterisk in Fig.~\ref{T1DIR}A and C). We assign this second kinetic signature to an intermediate state that directly precedes the formation of Pr.

\begin{figure*}[t]
	\centering
	\includegraphics[width=\linewidth]{T2DIR.png}
	\caption{Transient 1D~IR (top row) and transient 2D~IR spectra at time-points (50~\textmu\text{s}, 1~ms, and 11~ms) at which the population of certain intermediates is maximized, see Fig.~\ref{T1DIR}C. Transient 1D~IR spectra are taken from the data set of Fig.~\ref{T1DIR}A.}
	\label{FigTransient2DIR}
\end{figure*}

Fig.~\ref{FigTransient2DIR}D, E and F show a series of transient 2D~IR spectra at time points 50~\textmu\text{s}, 1~ms, and 11~ms, at which the populations of certain intermediates are maximized (compare dashed vertical lines in Fig. \ref{T1DIR}), together with the corresponding  transient 1D~IR spectra (Fig.~\ref{FigTransient2DIR} A, B, and C). Just like for transient 1D~IR spectroscopy (Fig.~\ref{2DIR}C, black vs. blue line), the late time transient 2D~IR spectrum at 11~ms is in essence the same as the stationary 2D~IR difference spectrum of Fig.~\ref{2DIR}F. All diagonal peaks of the latter can be identified also on the diagonal of the 11~ms transient 2D~IR  spectrum. In addition, cross peaks are revealed, but since the signal-to-noise of the 11~ms  transient 2D~IR  spectrum is not quite as good as that of the stationary 2D~IR difference spectrum, we focus the discussion here only on the strongest cross peak \textbf{III} between  CC$_\mathrm{AB}$ and CC$_\mathrm{BC}$. While the diagonal peak patterns are very similar in the earlier 1~ms transient 2D~IR  spectrum, that cross peak is not yet prominent. Hence, the coupling  between the corresponding two modes is switched on only during last kinetic step of the transition from Pg to Pr.

The 50~\textmu\text{s} transient 2D~IR  spectrum is very different. It has three diagonal peaks with inverted signs, in agreement with the corresponding bleaches in the transient 1D~IR spectrum at the same delay time (Fig.~\ref{FigTransient2DIR}A). Based on the results of the QM/MM calculations, we assign the two lower-frequency bleaches to the two bridging modes CC$_\mathrm{BC}$ and CC$_\mathrm{AB}$, while CC$_\mathrm{CD}$ is expected to be too weak to significantly contribute. Since the spectrum is mostly a bleach without any overlapping positive contribution, one  can read off the frequencies of the corresponding modes in the Pg state, listed in Table~\ref{Table}.
Furthermore, an inverted cross peak \textbf{II} can  be identified between CC$_\mathrm{BC}$ and CC$_\mathrm{AB}$. That is, as the bands on the diagonal bleach, the corresponding cross peaks bleaches as well. The calculations predict the strongest coupling between the adjacent  CC$_\mathrm{AB}$ and CC$_\mathrm{BC}$ modes (cross peak \textbf{II} in Fig. \ref{2DIR}J), in agreement with the present assignment of the  50~\textmu\text{s} transient   spectra.

Based on the QM/MM results, we also conclude that third negative band at around 1655~cm$^{-1}$ cannot originate from the chromophore, hence it must originate from the protein. Indeed, in our previous study,\cite{Buhrke2020b} we observed that the kinetic process leading to that band on a 40~ns timescales is only observed in the IR transient spectra, and not in the electronic spectra of the chromophore, which lead us to  conclude that it reflects a response of the protein. Cross peak \textbf{VI} indicates the coupling between the chromophore and the protein (Fig.~\ref{FigTransient2DIR}A).

\section{Discussion}

In transient IR spectroscopy of photoactive proteins, signals that originate from the chromophore are often dominant, since its electronic and geometric structure changes significantly more than that of the embedding protein.\cite{Gerwert1990, Ritter2015, Kottke2017, Lorenz-fonfria} Notable exceptions are phytochromes and LOV domains, where a large-scale re- or unfolding of entire secondary structure elements occurs.\cite{Stojkovic2014, Konold2016} The x-ray structures of Slr-g3 shown in Fig.~\ref{intro} reveal that the biggest changes in this protein concern rings A and D of the PCB chromophore. While all four rings are almost coplanar in the Pr state, rings A and D are tilted out of the plane in Pg. These geometries immediately explain the color change between the two states, i.e., their electronic properties: When rings A and D become coplanar with rings B and C, they increase the size of conjugated system of the chromophore, thereby shifting its electronic absorption bands towards the red.\cite{Wiebeler2018} Additionally, some differences in secondary structure have also been identified by x-ray crystallography, such as e.g. the unfolding of a small segment of the so-called a3-helix in the Pr state.\cite{Xu2020}

The kinetic analysis of our time-resolved data indicates two processes in the \textmu\text{s} to ms time range, or in other words, two meta-stable intermediate states that are maximized around 50~\textmu\text{s} and 1~ms, respectively, before Pr is formed. Both of these meta-stable configurations lie within the time resolution of the transient 2D~IR experiment, and thus the connectivity of the dominant modes allows unprecedented insight into the chromophore configuration at this stage of the photocycle.

The 50~\textmu\text{s} intermediate (Fig. \ref{FigTransient2DIR}A and B) only shows bleach bands associated with the conjugated system of PCB and almost completely lacks positive counterparts, indicating that the chemical structure of the PCB chromophore in this species differs significantly form all other intermediates and parent states. One possibility is that the PCB chromophore is extremely disordered, leading to very broad positive bands with weak intensity. Another process that could explain such spectra is transient deprotonation of the conjugated system, which leads to a reduced dipole strength and therefore reduced IR intensities of associated bands. Such a behaviour is well known from retinal model compounds as well as retinal-binding proteins, where protonated (positively charged) Schiff base compounds typically produce 5-10x stronger IR intensities in the C=C stretching region than their deprotonated (neutral) counter parts that have a lower level of charge alternation between the conjugated carbon atoms.\cite{Siebert1980, Braiman1987, Vogel2001} Transient de- and re-protonation events are also typical for the photocycle of phytochromes\cite{VanThor2001,Borucki2005}, and some studies identified such processes in the red/green CBCR AnPixJ.\cite{Escobar2013, Song2015} Interestingly, pH-dependent experiments on the same sample indicated that the protonation state of the chromophore has a large affect on the vibrational spectra, but barely affects the visible absorption properties\cite{Song2015, Altmayer2022}, which might explain why the  250~\textmu\text{s} process could not be resolved in the transient visible experiments.\cite{Xu2014} According to NMR and RR experiments, the PCB chromophore deprotonates either at the ring B or C pyrrole nitrogen, resulting in a tautomeric species where only one proton is shared between these two pyrrole rings.\cite{Song2015, Altmayer2022} In the light of these findings, we assign the 50~\textmu\text{s} intermediate to such a deprotonated PCB species.

The next intermediate is maximized at 1~ms, and shows positive counterparts to the bleach bands, confirming that the 250~\textmu\text{s} process indeed corresponds to PCB reprotonation. Here, the transient visible study suggested that the formation of Pr in the early milliseconds is preceded by an orange-absorbing intermediate state,\cite{Xu2014} and a RR study found that ring D is already in a Pr-like configuration while ring A still retains an out-of-plane tilt similar to Pg.\cite{Buhrke2020} Such a configuration has a conjugated system with an size that is in-between that of Pg and Pr, hence an electronic absorption band between green and red. The transient 1D~IR results of Fig. \ref{T1DIR} validate this hypothesis, because here it can be directly observed that the final difference signal associated with CO$_\mathrm{D}$ at 1674~cm$^{-1}$ is established at an earlier point in time than CO$_\mathrm{A}$ at $>$1700~cm$^{-1}$. Interestingly, The CO$_\mathrm{D}$ bleach signal is not accompanied by a visible positive counterpart until the reprotonation step, which means that the final coordination environment of ring D that defines the C=O frequency is only established during the 250~\textmu\text{s} process.

As ring A becomes coplanar during the last reaction step, the coupling between between CC$_\mathrm{AB}$ and CC$_\mathrm{BC}$ is switched on, as evidenced by cross peak \textbf{III} (Fig.~\ref{FigTransient2DIR}E vs F). No X-ray structures exist currently for any of the intermediates, and hence anharmonic QM/MM calculations as for Pg or Pr are not feasible. On a less rigorous level, it has been argued that the cross anharmonicity that determines the cross peak intensity in a 2D~IR spectrum correlates with the degree to which normal modes are delocalized, since the source of anharmonicity is local, i.e., that of individual chemical bonds.\cite{Woutersen2002,hamm11,Baiz2011} Upon a normal mode transformation, which delocalizes local modes in a harmonic sense, the local anharmonicity delocalizes accordingly. The complementary intensities of CC$_\mathrm{AB}$ and CC$_\mathrm{BC}$ in the Raman (Fig.~\ref{Raman}) vs the IR spectra (Table~\ref{Table}) indicate that the modes are delocalized to a certain extent in Pr. It it reasonable to assume that these modes are more local in a distorted structure with ring A tilted out of the plane, explaining the weaker anharmonic coupling.

\section{Conclusions}

A combination of three advanced IR methodologies allowed us to track the structural changes of the PCB chromophore inside the protein matrix of the red/green CBCR Slr-g3. First, transient 1D~IR spectroscopy was used to establish the time scales of various photocycle events. The higher resolution power of IR spectroscopy, as compared to UV/Vis spectroscopy, revealed an additional intermediate
preceding the formation of Pr, which has not been observed before. In addition, difference cross-peaks in the stationary 2D~IR spectrum guided the assignment of certain signals based on changes in vibrational coupling of the chromophore  modes. Here it is important to note that the discussed vibrational modes overlap with the strongly overwhelming amide I modes from the protein backbone. Nevertheless, the vibratonal modes of the chromophore can be isolated dy difference spectroscopy and validated by the calculations. Finally, the two previous experiments were combined into transient 2D~IR spectroscopy. The 2D~IR difference spectra of various intermediate species were obtained, that allowed us to gain insights into the coupling patterns associated with the later stages after the photoreaction. The combination of these three techniques allowed us to draw consistent conclusions about the reaction mechanism that Slr-g3 undergoes after excitation. For example, reaction steps could be determined during which the PCB chromophore de- and reprotonates and during which ring A rotates into plane to form the Pr state.

The work also demonstrates that a close-to-quantitative agreement between 2D-IR spectra and anharmonic QM/MM calculations  can be obtained for the initial and final forms, Pg and Pr, respectively. The existing X-ray structures of these two forms served as a basis of the QM/MM simulations. Here, we also make suggestions for the structurual changes of the two transient intermediates preceding the formation of final Pr state, based on a combination of  transient 1D and 2D spectroscopy.

Time-resolved X-ray scattering experiments are nowadays within reach, using upcoming X-ray free-electron laser facilities,\cite{Standfuss2019} which could provide the fundament of rigorous QM/MM calculations also of those intermediates, to test our structural hypotheses. Since Slr-g3 (as well as other CBCRs) are emerging tools in biotechnological applications, such as optogenetics and fluorescence microscopy,\cite{Oliinyk2019,Blain-Hartung2018} it is highly desirable to understand the molecular mechanism of its photoreaction in real time and at full atomistic detail.\\

\noindent \textbf{Acknowledgement:} This work was supported by the Swiss National Science Foundation (SNF) through Grant No. 200020B 188694/1. DB acknowledges a Liebig-Scholarship by the funds of the German chemical industry (Fonds der chemischen Industrie, FCI). I.S. acknowledges support from the Israel Ministry of Science and Technology (Grant 3-16311). I.S. thanks the DFG Collaborative Research Center 1078, project C6 for support.\\

\noindent \textbf{Data availability statement:} The data that support the findings of this study will be deposited on a repository at the proof stage.\\