# An acousto-optic modulator based bi-frequency interferometer for quantum
technology

Wenqi Li  College of Precision Instrument and Opto-Electronics Engineering,
Key Laboratory of Opto-Electronics Information Technology, Ministry of
Education, Tianjin University, Tianjin 300072, People’s Republic of China
Qiqi Deng  College of Precision Instrument and Opto-Electronics Engineering,
Key Laboratory of Opto-Electronics Information Technology, Ministry of
Education, Tianjin University, Tianjin 300072, People’s Republic of China
Xueshi Guo<EMAIL_ADDRESS>College of Precision Instrument and Opto-
Electronics Engineering, Key Laboratory of Opto-Electronics Information
Technology, Ministry of Education, Tianjin University, Tianjin 300072,
People’s Republic of China
Xiaoying Li<EMAIL_ADDRESS>College of Precision Instrument and Opto-
Electronics Engineering, Key Laboratory of Opto-Electronics Information
Technology, Ministry of Education, Tianjin University, Tianjin 300072,
People’s Republic of China

###### Abstract

We demonstrate a high performance AOM based bi-frequency interferometer, which
can realize either beating or beating free interference for single photon
level quantum state. Visibility and optical efficiency of the interferometer
are $(99.5\pm 0.2)\%$ and $(95\pm 1)\%$, respectively. The phase of the
interferometer is actively stabilized by using dithering phase locking scheme,
where the phase dithering is realized by directly driving the AOMs with
specially designed electronic signal. We further demonstrate applications of
the interferometer in quantum technology, including bi-frequency coherent
combination, frequency tuning and optical switching. These result show the
interferometer is a versatile device for multiple quantum technologies.

##  I Introduction

Making use of acousto-optic bragg diffraction effect, Acousto-optic modulators
(AOMs) are versatile active optical devices that can efficiently change both
the frequency and the propagating direction of an optical field in real time
[1]. AOMs have been used in various applications, such as stabilizing the
output power, reducing the line-width and shaping the spatial profile of a
laser beam etc. [2, 3, 4]. AOMs can be exploited to realize active optical
tweezers and ultra short pulse picking [5, 6]. In addition to the
aforementioned applications in classical technology, AOMs are also
indispensable for a lot of quantum technologies. Taking the advantage of
switching with high isolation, they have been used in photon subtraction based
Non-Gaussian state generation [7, 8], photon triggered homodyne tomography [9]
or controlling of a quantum memory [10]. Taking the advantage of introducing
frequency shifting to an optical field, they are used to implement optical
heterodyning [11, 12], to observe a beating signal from single photons [13] or
to generate a phase locking reference without displacing the quantum state
[14]. Currently, the maximum diffraction efficiency of commercially available
AOMs is around 85$\%$, which reduces as the bandwidth of AOM increases, and
harsh requirement on the input beam size and driving power has to been
fulfilled in order to achieve this efficiency. While the efficiency is enough
for most classical applications, it is still not high enough for the quantum
applications where quantum states have to transmit through the AOMs and keep
the high purity. This impels us to consider the interferometric enhancement of
the diffraction effect when using AOMs in quantum technology.

In most application schemes of AOMs, only one output port (either the direct
pass beam or the diffracted beam) is used and the other is abandoned [15, 16].
However, AOMs are apparently 4-port devices that can be used for bi-frequency
beam splitting. Recently, AOM is used as a beam-splitter to implement photon
heterodyning [12]. More recently, AOM-based interferometers are used to
observe the beating signal of single photons [13] or to bidirectional routing
laser pulses [17], where both beam-splitting and beam-combining functions are
carried out by AOMs. The visibility of the beating signal in Ref. [13] is only
about $15\%$, far from the ideal value of 1. Moreover, both experiments make
use of co-propagation optical configuration and create either slow beating of
about 1 kHz [13] or fast beating signal of about 200 MHz [17]. On one hand,
this co-propagation setup and scanning effect of beating signal will reduce
influence the phase drifting and make active phase stabilizing less important.
On the other hand, beating free and spatially separated interference can be
necessary in many application schemes. Here we demonstrate an AOM based bi-
frequency interferometer (ABI) working in chopped phase locking [18] mode,
which can realize either beating or beating free interference for single
photon level quantum state. By carefully designing the optical system and
optimizing the mode matching, the visibility of the beating signal is
$(99.5\pm 0.2)\%$. We also carefully optimize the optical transmission
efficiency of the system. Especially, with special designed driving signal for
the AOM, the modulation for phase locking [19] is carried out with the same
AOM for beam splitting. This simplifies the optical setup and further enhances
the efficiency of the system. Moreover, we propose the applications of the
interferometer in quantum technology, including bi-frequency coherent
combination, frequency tuning and optical switching. With the support of
experimental testing, we show these applications schemes can in principle
achieve near perfect efficiency while only about $50\%$ diffraction efficiency
is needed for the AOMs. This property will merits many applications for it
will greatly reduce the demand for the AOM driving and therefore enhance the
bandwidth of the system. These characteristics make the interferometer a
powerful tool for many quantum technologies.

Figure 1:  (a) The optical frequency shifting and interference property when
an acousto-optic modulator (AOM) is used as a bi-frequency beam
splitter/combiner. (b) An AOM based bi-frequency interferometer (ABI) whose
output ports preserve the optical frequency at the output ports

##  II Working principle of the ABI

We first illustrate the working principle of the ABI. As shown in Fig. 1(a),
without loss of the generality, we assume two optical fields
$\hat{a}_{\omega_{1}}$ and $\hat{b}_{\omega_{2}}$ are aligned to the -1st and
the +1st order Bragg diffraction of an AOM, and the output fields $\hat{c}$
and $\hat{d}$ have the form of [13]

$\displaystyle{\hat{c}}$ $\displaystyle=$ $\displaystyle t\hat{b}_{\omega
2}+e^{i\theta}r\hat{a}_{\omega 1-\Omega}$ $\displaystyle{\hat{d}}$
$\displaystyle=$ $\displaystyle t\hat{a}_{\omega
1}-e^{-i\theta}r\hat{b}_{\omega 2+\Omega},$ (1)

where $\hat{a}$ and $\hat{b}$ denote the spatial modes and $\omega_{1(2)}$ in
the subscripts denotes their optical frequency. After passing the AOM, the
optical frequency is up/down shifted by the frequency of the radio frequency
(RF) signal $\Omega$. $t$, $r$ satisfying $t^{2}+r^{2}=1$ are the real number
transmission and reflection parameters determined by the strength of the RF
signal, and $r^{2}$ is the diffraction efficiency of the AOM. $\theta$ is AOM
induced phase decided by the phase of the RF signal. Each field operator in
Eq. (II) satisfies the standard bosonic commutation relation, for example
$\hat{a}$ has the commutation relation $[\hat{a},\hat{a}^{\dagger}]=1$.

Assuming the frequency of ports $\hat{a}$ and $\hat{b}$ are $\omega+\Omega$
and $\omega$, the ouput ports $\hat{c}$ and $\hat{d}$ only contains the
optical field with frequency of $\omega$ and $\omega+\Omega$, respectively.
Therefore, $\hat{a}_{\omega+\Omega}$ and $\hat{b}_{\omega}$ interfere at
$\hat{c}_{\omega}$ and $\hat{d}_{\omega+\Omega}$, and can be coherently
combined by using another AOM to form a bi-frequency Mach-Zehnder
interferometer as shown in Fig. 1(b). The output fields of the interferometer
can be obtained by using Eq. (II) twice:

$\displaystyle{\hat{e}}_{\omega+\Omega}=t_{1}^{\prime}{\hat{a}}_{\omega+\Omega}-r_{1}^{\prime}{\hat{b}}_{\omega+\Omega}$
$\displaystyle{\hat{f}}_{\omega}=r_{2}^{\prime}{\hat{a}}_{\omega}+t_{2}^{\prime}{\hat{b}}_{\omega}$
(2)

with

$\displaystyle
t^{\prime}_{1}=t_{1}t_{2}e^{i\varphi}-r_{1}r_{2}e^{i(\theta_{1}-\theta_{2})}$
$\displaystyle
r^{\prime}_{1}=r_{1}t_{2}e^{i(\varphi-\theta_{1})}+t_{1}r_{2}e^{-i\theta_{2}}$
$\displaystyle
t^{\prime}_{2}=t_{1}t_{2}-r_{1}r_{2}e^{i(\theta_{2}-\theta_{1}+\varphi)}$
$\displaystyle
r^{\prime}_{2}=r_{1}t_{2}e^{i\theta_{1}}+t_{1}r_{2}e^{i(\theta_{2}+\varphi)}$
(3)

These effective coefficients are decided by the optical path induced phase
difference $\varphi$, the AOM splitting parameters $t_{1(2)}$, $r_{1(2)}$, and
the AOM driving induced phase $\theta_{1(2)}$ for AOM1 and AOM2, respectively.

Assuming we send light from $\hat{b}$ port and detect it from $\hat{e}$ port
in Fig. 1(b). Then when the input light is a coherent state $|\alpha\rangle$
and t, r are set to $t_{1(2)}^{2}=r_{1(2)}^{2}=0.5$, the intensity at
$\hat{e}$ can be obtained by using Eq. (II)

$\mathrm{I}_{out}=\langle\alpha|{\hat{e}_{\omega+\Omega}^{\dagger}}{\hat{e}_{\omega+\Omega}}|\alpha\rangle=\frac{1}{2}\Big{[}1+\cos({\phi})\Big{]}\mathrm{I}_{in},$
(4)

where $\phi=\varphi-\theta_{1}+\theta_{2}$ is the overall phase parameter of
the interferometer, $\mathrm{I}_{in}=|\alpha|^{2}$ is the intensity of the
input light. Eq. (4) describes the ABI in ideal condition. In practice, there
inevitably exists the three imperfections: (i) the non-ideal efficiency of the
optical components, (ii) the non-ideal mode matching between the two arms of
the interferometer , and (iii) the frequency difference between the two RF
driving signals. With these factors considered, Eq. (4) is rewritten as

$\mathrm{I}_{out}^{\prime}=\frac{\eta}{2}\Big{[}1+V\cos({\Delta\omega
T+\phi})\Big{]}\mathrm{I}_{in},$ (5)

where $\eta$ is the efficiency determined by the losses of the optical
elements in ABI, $V$ is the mode mis-matching induced visibility of the
interference, and $\Delta\omega T$ is the beating term induced by the
frequency difference of the RF signals that drive the two AOMs.

##  III Applications in quantum technology

To show the ABI with high efficiency and high visibility is a very useful tool
for quantum technology, we analyze three typical functions that can be
implemented by using ABI (see Fig. 2). As shown in Fig. 2(a), the ABI can be
applied to coherently combine two quantum state with a slightly different
frequencies. From Eq. (II-II), one sees the ABI can achieve arbitrary
effective splitting ratio $t^{\prime}_{1(2)}$ and $r^{\prime}_{1(2)}$ when the
maximum diffraction efficiency of each AOM can reach $50\%$ and the phase
parameters $\varphi$, $\theta_{1}$ and $\theta_{2}$ are well controlled. Since
frequency is one of the important degree of freedoms for optical mode [20],
this scheme can be used to create continuous variable Einstein-Podolsky-Rosen
state out of two single mode squeezing state with frequency difference or
creating a large scale entangled state with multiple pairs of frequency
encoded two mode entangled states [21, 22, 23].

Figure 2:  Three typical functions fulfilled by the AOM based bi-frequency
interferometer (ABI). (a) A Bi-frequency beam splitter with variable splitting
ratio. (b) A frequency tuner of quantum state. (c) An optical switch
controlled by gating the RF driving signal to have the AOMs chopped

The scheme in Fig. 2(a) can interfere two quantum states with different
frequency. By slightly modifying the scheme (see Fig. 2(b)), the ABI can
realize a frequency tuner scheme for quantum states, which can tune the
frequency of a quantum optical state and fit that of the other quantum system
[24, 10]. In this scheme, a quantum state $|\psi\rangle_{\omega}$ is sent in
at $\hat{b}$ port and the $\hat{a}$ port is reserved for vacuum input.
According to Eq. (1-4), by controlling the overall phase parameter of the ABI
to $\phi=0$, the quantum state will be shifted by $\Omega$ in optical
frequency and output at $\hat{e}$ port, where $|\psi\rangle_{\omega+\Omega}$
is used to denote the output state with frequency shifted. In this ideal case,
the state at $\hat{f}$ port is a vacuum state. Since this frequency shift is
controlled by the frequency $\Omega$ of the RF signal, it provide us a handy
tool to realize quantum state frequency tuning.

If we control the overall phase parameter of the ABI to $\phi=0$ and $\pi$ in
Fig. 2(b), the input state $|\psi\rangle_{\omega}$ will go to port $\hat{e}$
with frequency shift and port $\hat{f}$ without frequency shift, respectively.
Therefore, the ABI can also be used as an optical switch. However, in this
scheme, both the isolation and the transmission efficiency of the switch is
confined by the visibility of the interference. In quantum technologies, AOM
based switches are often used as a protection stage before a week light
detection process [7, 9], and good isolation is a critical requirement for
such a switch. To solve this problem, we propose the optical switch scheme in
Fig. 2(c), where high isolation property is achieved by adding a gate to the
RF signals that drive the AOMs. When the RF signals are off, the optical path
of the interferometer is degraded into two independent propagated spatial
modes. With the same input configuration as Fig. 2(b), the input state will go
to $\hat{f}$ port when RF signal is off and the isolation between $\hat{e}$
and $\hat{f}$ port is high. When the RF signal is on, we are able to set the
effective transmission/reflection coefficients to $t^{\prime}_{1}=0$ and
$r^{\prime}_{1}=1$ by controlling the overall phase $\phi=0$. In this case,
the quantum optical state is transmitted to $\hat{e}$ port with frequency up
shifted.

We note all three schemes in Fig. 2 can also be implement with a single AOM.
However, benefiting from the interferometric enhancement of the diffraction
effect, our scheme can achieve a near perfect efficiency when each AOM is
working at around 50$\%$ diffraction efficiency. This enhancement is critical
in many quantum application schemes, especially for those using continuous
variable quantum state.

##  IV Experimental demonstration of the ABI

Our experimental setup of the acousto-optic modulator based bi-frequency
interferometer is shown in Fig. 3(a). A linearly polarized continuous wave
(CW) laser beam with a wavelength of 1550 nm is coupled to the +1st mode
($\hat{b}_{\omega}$ in Fig. 1(b)) of the AOMs (ISOMET M1205) with a fiber
collimator (FC). The beam is split into a frequency up-shifted arm (red beam
with the frequency of $\omega+\Omega$) and a non-shifted arm (blue beam with
the frequency of $\omega$) by AOM1. The two beams are then coherently combined
by AOM2. For each arm between two AOMs, two plane mirrors and a concave mirror
with 300 mm radius of curvature are used to ensure the required mode matching
is optimized. To control the phase between two arms, a piezo device (PZT) is
mounted on the concave mirror of the frequency shifted arm. At the output of
the AOM2, the frequency up-shifted port is detected by power detector, which
is either a photon diode (PD1) for measuring intensity with higher photon
number or a single photon detector (SPD, id Quantique-id200) for intensity at
single photon level. The output of the frequency non-shifted port is detected
with PD2. Detection result for both PD1/SPD and PD2 are sent to a FPGA
(STEMlab125-14) based data-acquisition system (DAQ), which also generates the
feedback signal to the PZT for phase locking.

Figure 3:  (a) The experimental setup of the AOM based bi-frequency
interferometer. AOM, Acousto-optic modulator; CM, concave mirror; PZT, piezo-
electric translation device; PD, photo diode; SPD, single photon detector;
DAQ, data acquisition system. (b) 100 kHz beating signal detected by PD1 from
the light with micro-Watt level intensity. (c) 5 Hz Beating signal detected
with SPD from the light with single photon level intensity. Red dash-dotted
lines in (b) and (c) are the fitting to the data  Figure 4:  (a) Output
intensity of an ABI without beating detected by PD1 for PZT scanned (blue),
locked to maximum (green) and lock to minimum (orange). (b) The AOM driving
signal RF1 and RF2 to realize phase dithering for the phase locking scheme.
Both RF1 and RF2 are centered at 80 MHz. The modulation for phase locking of
the ABI is realized by adding a 200 kHz small signal phase modulation to RF1.
For the illustration propose, the upper trace is plotted with larger
modulation strength and frequency than it is used in the experiment

We first run the setup in slow beating mode, where the active phase
stabilizing is off and no signal is added on the PZT. In this mode, we
demonstrate the bi-frequency property by the beating effect, and use the beat
result to calibrate the visibility $V$ and the transmitting efficiency $\eta$
of the system. RF1 and RF2 drives the AOMs are generated by a two channel
arbitrary function generator (Siglent-SDG2122X), and are set to 80 MHz and
79.9 MHz. Thus the beating frequency $\Delta\omega/2\pi$ is 100 kHz. Fig. 3(b)
shows a measured beating signal for the input light of micro-Watt level
intensity, where $\mathrm{I}^{\prime}_{out}$ is normalized to the intensity of
input field $\mathrm{I}_{in}$. The blue trace are the directly measured data
and the red dash-dotted line are the fitting to Eq. (5). $V$, $\eta$ and
$\phi$ in Eq. (5) are set to be the fitting parameters. One sees the fitting
curve overlaps well with the directly measured data, and the fitting result is
V=$(99.5\pm 0.2)\%$, $\eta$=$(95\pm 1)\%$ and $\phi$= (1.08$\pm$0.01) rad,
respectively. This result shows we realize almost perfect visibility and high
optical efficiency simultaneously. We note each optical surface in our setup
has a loss of about $0.5\%$ limited by the quality of optical coating, so
$\eta$ can be further improved once components with better optical coatings
are used. We then heavily attenuate the intensity of the input light to single
photon level, and detect the beating signal by replacing PD1 with the SPD,
which is triggered by its internal 50MHz clock. To get enough counts for a
counting window (10 ms), the beating frequency is reduced to 5 Hz by changing
the driving frequency of AOM2. The counting result and its fitting are shown
in Fig. 3(c), where the detection efficiency of the SPD is about $20\%$ and
the dark count is $(4\times 10^{-6})$ counts/trigger at this counting setup.
Compared with the maximum counting rate of about 0.06 counts/trigger, the dark
count is negligible in this measurement. Compared to the result in Fig. 3(b),
which is obtained in the condition of input light at mico-watt level and
beating frequency of 100 KHz, data points in Fig. 3(c) obviously deviate from
the fitting curve due to larger time scale for phase drifting. However, this
will not prevent us to get the visibility since it is only determined by the
maximum and the minimum value of the curve. The fitting result shows the
visibility for light with optical intensity of single photon level is
$(99.2\pm 0.3)\%$, which is almost the same as the case for input light with
higher intensity.

Secondly, we set RF1 and RF2 to be identical frequency of exactly 80 MHz and
scan the PZT with a 30 Hz ramp signal. The interference fringe detected by PD1
is shown by the blue trace in Fig. 4(a), where the input light is on the
micro-Watt level. When triggered by the ramp signal, one sees the fringe moves
on the oscilloscope horizontally in random direction. This is caused by the
phase drift between the two arms, and this show it is necessary to introduce
active phase locking scheme. We use dithering phase locking scheme [19] to
stably lock the ABI to arbitrary phase. As it is shown in Fig. 4(b),
modulation required in this phase locking scheme is realized by adding a 200
kHz phase modulation with a small modulation depth to RF1. Therefore, AOM1 is
used for both beam splitting and introducing of the modulation for phase
locking. This special design avoids using extra optical modulators and
enhances the optical efficiency of the system. By properly configuring the
FPGA system to implement digital demodulators, PID controllers and digital
filters [25], the overall phase $\phi$ can be locked to arbitrary value using
the detection result of PD1(2). The green and orange traces in Fig. 4(a) shows
two typical locking results when the overall phase $\phi$ is $0$ and $\pi$ and
the values of $\mathrm{I}_{out}/\mathrm{I}_{in}$ at $\hat{e}$ port is locked
to the maximum and minimum point, respectively. Compare the beating signal, we
find the visibility of interference for the identical frequency case is
reduced to $(93.7\pm 0.5)\%$. Substituting the $\eta$ obtained from the
fitting of the beating signal and this visibility into Eq. (5), we find, the
total transmission efficiency of the ABI when used as an optical switch is
$\mathrm{I^{\prime}_{out}/I_{in}=}(92\pm 1)\%$. We study the reason for this
visibility reduction and find, besides a phase shift, the PZT can cause the
beam walking-off and induce mode mis-matching. This visibility reduction can
be in principle avoid if the feedback signal is applied to RF2, and AOM2 is
directly used as the phase shifting device of our phase locking system instead
of the PZT. To complete this optimization, angle calculation and phase
unwrapping algorithms is necessary to be implemented in the FPGA [24] and the
progress is underway.

Figure 5:  (a) Output intensity in switch mode detected by PD1 with a locking
duty cycle of $30\%$. In this mode, the AOM driving signal RF1(2) are turned
on and off periodically by a gate signal shown above the intensity result. (b)
SPD counting result of the ABI in frequency tuner mode with a detection duty
cycle of $50\%$. The error bar for counting result is the same size or smaller
than the data point. The input light is a temporally complementary pattern of
a locking reference (LR) light with intensity on the micro-Watt level and a
week coherent state (Coh) with intensity on the single photon level. The
timing sequence of LR and Coh are shown above the counting result. SPDe, SPD
trigger enable signal

In the next, we run the setup in a chopped locking mode for a gated RF signal.
In this experiment, the input light is CW with intensity on the micro-Watt
level, but the driving signal RF1 and RF2 are gated by two voltage controlled
attenuators. This operating mode realizes the optical switch scheme in Fig.
2(c). The RF gate signal is 100 Hz in repetition rate and the duty cycle is
tunable. When the gate signals are on high voltage, the whole system works the
same as ABI being locked at the maximum output. When the gate signals are low,
the driving signal of the two AOMs in the ABI is completely cut off, and the
injected light will go through the frequency non-shifted arm (blue beam) in
Fig. 3(a), leaving the frequency shifted arm (red beam) in high isolation
state to the input field. A locking feedback enable signal synchronized to the
RF gate signal (see the yellow waveform labeled by ’Gate to RF1(2)’ in Fig.
5(a)) is send to the DAQ to stop the locking process and hold the feedback
value to increase the stability of the locking process. In the experiment we
gradually reduce the duty cycle of the RG signal and find the setup works well
with a duty cycle as low as 30$\%$, and a typical measured intensity result
detected by PD1 is shown by the blue trace in Fig. 5(a), in which one sees the
lock are stable after a locking establishing time of sub-milli-second order.
We note the phase locking stability increases when either the frequency or the
duty cycle of the locking process increases. We further measure the isolation
of the ABI when used as an optical switch. Since the isolation exceeds the
dynamic range of the SPD, we insert calibrated optical attenuators after the
FC when measuring the photon numbers for RF on state and remove it when
measuring RF off state. By comparing the photon counting result for the on and
off state, we find the isolation is about $74$ dB. This shows the setup
functions as a high isolation optical switch from either input port to its
double diffracted output port. Though confined by the visibility, isolation of
either input port to its double direct pass port is only around 15 dB when
used as an optical switch.

Finally, we run the setup in a chopped phase locking mode for a temporally
complementary input of a locking reference (LR) light and a week coherent
light (Coh). The high extinction ratio of $\geq 110$ dB between LR and Coh is
ensured by using two mechanical fiber optical switches (e-Photics YFS-2x2),
which prevent the SPD from counting the photons in LR. This mode realizes the
frequency tuner scheme in Fig. 2(b) for a week input state. The waveform at
the top of Fig. 5(b) shows the timing sequence of this mode, where the micro-
Watt level LR is on for $30\%$ of the time and has a frequency of 5 Hz. LR is
detected at the frequency non-shifted output port by PD2 and generate feedback
signal for phase locking, which is also enabled by a locking feedback enable
signal synchronized to LR. When LR is off, the feedback signal is hold and the
RF signals continue to drive the AOMs. Coh with a intensity on the single
photon level is used as input to simulate a quantum state
$|\psi\rangle_{\omega}$ and detected at the frequency up-shifted port by the
SPD. To avoid saturating the SPD, we only enable the SPD using a gate with
$50\%$ duty cycle by using a SPD enable signal SPDe. Data points in Fig. 5(b)
shows the counting result with a counting window of 20 ms, where one sees the
counting rate stability follows the shape of SPDe signal. This result shows
our ABI can work in a chopped locking scheme for single photon level input.

##  V Summary

In summary, we experimentally realize an AOM based bi-frequency
interferometer. The interference visibility of the beating signal of the ABI
is (99.5$\pm 0.2)\%$ and the optical efficiency is $(95\pm 1)\%$. Benefiting
from the interferometric enhancement of the diffraction effect, the
interferometer can realize multiple quantum technologies in high efficiency
such as bi-frequency coherent combining, optical switching, and frequency
tuning. The interferometer can work in chopped locking mode, which enables the
operating of light with intensity on the single photon level. With specially
designed RF driving signal, the modulation for phase locking is implemented on
the beam splitting AOM, which reduce the number of optical elements and
increase the overall quantum efficiency. The performance of experimental
system we realized is currently confined by the optical loss of the components
and the beam walk-off induced by the PZT. With the high visibility of the
beating signal, we believe the overall efficiency to implement optical
switching or coherently combining of quantum states with the ABI can be
further improved to around $98\%$-$99\%$ by using commercially available
ultra-low loss optical components.

###### Acknowledgements.

This work was supported in part by National Natural Science Foundation of
China (Grants No.12004279 and 12074283).

## References

  * Wade [1990] G. Wade, Bulk-wave acousto-optic bragg diffraction, in  _Guided-Wave Acousto-Optics: Interactions, Devices, and Applications_ , edited by C. S. Tsai (Springer Berlin Heidelberg, Berlin, Heidelberg, 1990) pp. 11–52.
  * Tricot  _et al._ [2018] F. Tricot, D. H. Phung, M. Lours, S. Guérandel, and E. de Clercq, Power stabilization of a diode laser with an acousto-optic modulator, Review of Scientific Instruments 89, 113112 (2018).
  * Dong  _et al._ [2015] J. Dong, Y. Hu, J. Huang, M. Ye, Q. Qu, T. Li, and L. Liu, Subhertz linewidth laser by locking to a fiber delay line, Appl. Opt. 54, 1152 (2015).
  * Song  _et al._ [2020] E. Song, T. Dai, G. Zhu, H. Wang, K. Aleksei, and X. Zhu, Adjustable and stable beam profile generation in a yb:yag thin-disk laser, Opt. Lett. 45, 6550 (2020).
  * Bola  _et al._ [2020] R. Bola, D. Treptow, A. Marzoa, M. Montes-Usategui, and E. Martín-Badosa, Acousto-holographic optical tweezers, Opt. Lett. 45, 2938 (2020).
  * de Vries  _et al._ [2015] O. de Vries, T. Saule, M. Plötner, F. Lücking, T. Eidam, A. Hoffmann, A. Klenke, S. Hädrich, J. Limpert, S. Holzberger, T. Schreiber, R. Eberhardt, I. Pupeza, and A. Tünnermann, Acousto-optic pulse picking scheme with carrier-frequency-to-pulse-repetition-rate synchronization, Opt. Express 23, 19586 (2015).
  * Takase  _et al._ [2022] K. Takase, A. Kawasaki, B. K. Jeong, M. Endo, T. Kashiwazaki, T. Kazama, K. Enbutsu, K. Watanabe, T. Umeki, S. Miki, H. Terai, M. Yabuno, F. China, W. Asavanant, J. ichi Yoshikawa, and A. Furusawa, Generation of schrödinger cat states with wigner negativity using a continuous-wave low-loss waveguide optical parametric amplifier, Opt. Express 30, 14161 (2022).
  * Liu  _et al._ [2022a] S. Liu, D. Han, N. Wang, Y. Xiang, F. Sun, M. Wang, Z. Qin, Q. Gong, X. Su, and Q. He, Experimental demonstration of remotely creating wigner negativity via quantum steering, Phys. Rev. Lett. 128, 200401 (2022a).
  * Kawasaki  _et al._ [2022] A. Kawasaki, K. Takase, T. Nomura, S. Miki, H. Terai, M. Yabuno, F. China, W. Asavanant, M. Endo, J. ichi Yoshikawa, and A. Furusawa, Generation of highly pure single-photon state at telecommunication wavelength, Opt. Express 30, 24831 (2022).
  * Appel  _et al._ [2008] J. Appel, E. Figueroa, D. Korystov, M. Lobino, and A. I. Lvovsky, Quantum memory for squeezed light, Phys. Rev. Lett. 100, 093602 (2008).
  * Michaud-Belleau  _et al._ [2022] V. Michaud-Belleau, J.-D. Deschênes, and J. Genest, Reaching the true shot-noise-limited phase sensitivity in self-heterodyne interferometry, IEEE Journal of Quantum Electronics 58, 1 (2022).
  * Okawa  _et al._ [2017] Y. Okawa, F. Omura, Y. Yasutake, and S. Fukatsu, Photon heterodyning, Opt. Express 25, 20156 (2017).
  * Mathevet  _et al._ [2020] R. Mathevet, B. Chalopin, and S. Massenot, Single photon beat note in an acousto-optic modulator-based interferometer, American Journal of Physics 88, 313 (2020).
  * Arnbak  _et al._ [2019] J. Arnbak, C. S. Jacobsen, R. B. Andrade, X. Guo, J. S. Neergaard-Nielsen, U. L. Andersen, and T. Gehring, Compact, low-threshold squeezed light source, Opt. Express 27, 37877 (2019).
  * Ma  _et al._ [2020] X. Ma, X. Zhang, K. Huang, and X. Lu, Noise-suppressing and lock-free optical interferometer for cold atom experiments, Opt. Express 28, 28584 (2020).
  * de Lima  _et al._ [2006] M. M. de Lima, M. Beck, R. Hey, and P. V. Santos, Compact mach-zehnder acousto-optic modulator, Applied Physics Letters 89, 121104 (2006).
  * Liu  _et al._ [2022b] R. Liu, Y. Ma, L. Ji, L. Qiu, M. Ji, Z. Tao, and S. Wu, Composite acousto-optical modulation, Opt. Express 30, 27780 (2022b).
  * Herzog  _et al._ [2006] F. Herzog, K. Kudielka, D. Erni, and W. Bachtold, Optical phase locking by local oscillator phase dithering, IEEE Journal of Quantum Electronics 42, 973 (2006).
  * Wu  _et al._ [2019] S. Wu, W. Huang, P. Yang, S. Liu, and L. Chen, Arbitrary phase-locking in mach–zehnder interferometer, Optics Communications 442, 148 (2019).
  * Chen  _et al._ [2014] M. Chen, N. C. Menicucci, and O. Pfister, Experimental realization of multipartite entanglement of 60 modes of a quantum optical frequency comb, Phys. Rev. Lett. 112, 120505 (2014).
  * Larsen  _et al._ [2019] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Fiber-coupled epr-state generation using a single temporally multiplexed squeezed light source, npj Quantum Information 5, 46 (2019).
  * Yoshikawa  _et al._ [2016] J.-i. Yoshikawa, S. Yokoyama, T. Kaji, C. Sornphiphatphong, Y. Shiozawa, K. Makino, and A. Furusawa, Invited article: Generation of one-million-mode continuous-variable cluster state by unlimited time-domain multiplexing, APL Photonics 1, 060801 (2016).
  * Asavanant  _et al._ [2019] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J. ichi Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Furusawa, Generation of time-domain-multiplexed two-dimensional cluster state, Science 366, 373 (2019).
  * Suleiman  _et al._ [2022] I. Suleiman, J. A. H. Nielsen, X. Guo, N. Jain, J. Neergaard-Nielsen, T. Gehring, and U. L. Andersen, 40 km fiber transmission of squeezed light measured with a real local oscillator, Quantum Science and Technology 7, 045003 (2022).
  * Neuhaus  _et al._ [2017] L. Neuhaus, R. Metzdorff, S. Chua, T. Jacqmin, T. Briant, A. Heidmann, P.-F. Cohadon, and S. Deléglise, Pyrpl (python red pitaya lockbox) - an open-source software package for fpga-controlled quantum optics experiments, 2017 European Conference on Lasers and Electro-Optics and European Quantum Electronics Conference , EA$\\_$P$\\_$8 (2017).


# Simulating short- and long-term evolutionary dynamics on rugged landscapes

Leonardo Trujillo1, Paul Banse1 Guillaume Beslon1

1Inria Beagle Team, Artificial Evolution and Computational Biology,
Université de Lyon, Inria, ECL, INSA Lyon, Université Claude Bernard Lyon 1,
Univ. Lumière Lyon 2,
CNRS, LIRIS UMR 5205, F-69622, France
{leonardo.trujillo, paul.banse<EMAIL_ADDRESS>
###### Abstract

We propose a minimal model to simulate long waiting times followed by
evolutionary bursts on rugged landscapes. It combines point and inversions-
like mutations as sources of genetic variation. The inversions are intended to
simulate one of the main chromosomal rearrangements. Using the well-known
family of NK fitness landscapes, we simulate random adaptive walks, i.e.
successive mutational events constrained to incremental fitness selection. We
report the emergence of different time scales: a short-term dynamics mainly
driven by point mutations, followed by a long-term (stasis-like) waiting
period until a new mutation arises. This new mutation is an inversion which
can trigger a burst of successive point mutations, and then drives the system
to new short-term increasing-fitness period. We analyse the effect of genes
epistatic interactions on the evolutionary time scales. We suggest that the
present model mimics the process of evolutionary innovation and punctuated
equilibrium.

## Introduction

A very simple model for artificial molecular evolution consists of genetic-
like information stored and encoded in a binary sequence. If this digital
genome is subject to variations caused by mutations, along with selection,
then a minimalist evolutionary dynamics can be simulated. Can such a simple
model recreate different evolutionary time scales in a single framework? This
is a question that is not devoid of interest to evolutionary biologists, that
has caused controversy when trying to understand how the causal relationship
between micro and macro evolutionary processes could be (Uyeda et al.,, 2011).
In this sense, in contrast to Darwinian gradualism, c.f. (Darwin,, 1964, p.
32), the theory of punctuated equilibrium (PE) proposes that macroevolutionary
dynamics consists of short periods of rapid evolution alternated by long
periods of stasis with little evolutionary change, e.g. (Gould and Eldredge,,
1977). However, some of the statistical methods used in the quantitative
analysis of the historical data of macroevolutionary patterns, and their
relationship with PE, have also been the subject of controversy and confusion
(Pennell et al.,, 2014). On the other hand, several theoretical and
computational models have addressed many features of PE. For example, two
recent attempts to link the micro and the macro processes are: the
developmental gene hypothesis proposed by Casanova and Konkel, (2020), that
supports the genetic basis of PE; and the threshold response to increasing
ecological pressures resulting from evolutionary diversification at an
existing level of complexity, leading to punctuated expansions (Ispolatov et
al.,, 2019).

Far from providing an ultimate answer to biologists, in this paper we report
simulations that can simultaneously take into account changes across a range
of well-differentiated time scales, without incorporating exogenous elements
beyond the molecular evolution of our digital genome model (such as
environmental contingencies). We use the Kauffman, (1993) NK model of
molecular fitness landscapes, where $N$ is the length of the digital genome
and $K$ is the “epistatic” coupling between the $N$-sites in the genome and
determines the fitness landscape ruggedness (Kauffman and Levin,, 1987;
Kauffman and Weinberger,, 1989; Kauffman,, 1993).

In this framework, we elucidate that the differences between the time scales
are related to the occurrence of chromosomal rearrangements (Griffiths et
al.,, 2012, Ch. 17) (at least those related to mutations that reverse the
order of the genomic sequence, as reported in this work). Moreover, from the
simulations outcomes we identify evolutionary bursts. That is, once a new
inversion promotes an increase in fitness – after a long waiting time – it
triggers a successive occurrence of point mutations. We show that in an early
stage of the evolutionary process the system had exhausted the increasing-
fitness point mutations around a local maximum of fitness. Then, once an
inversion accesses a “new” domain of the fitness landscape, the exploration is
re-started and is carried out mostly by point mutations. Here, the meaning of
“new” refers to domains of the fitness landscape that would be inaccessible
through point mutations once a local maximum fitness value has been reached.
This provides a very simple and evocative notion to think about evolutionary
innovation and punctuated equilibrium.

In the following section of this paper we introduce the model, with a
description of our double-stranded representation of binary genomes. Then we
derive a new inversion-like mutation operation, together with the standard
bit-flip mutations. We continue with a brief introduction of the NK model and
the adaptive walk approach. Next, we present our simulations and our main
results. We close the paper with some conclusions suggesting that the present
model mimics the process of evolutionary innovation and paves the way to
simple simulations to explore dynamics such as punctuated equilibrium.

## The model

### Genome representation

We assume that the digital organism’s genotype is a Boolean sequence of length
$N\in\mathbb{N}$. In this way, all possible genome combinations define a
$2^{N}$ Boolean hypercube space $\mathcal{Q}_{N}=\\{0,1\\}^{N}$. Inspired by
molecular biology, these genotypes are double-stranded binary vectors ${\bf
x}:=(x_{1},\ldots,x_{N})\in\mathcal{Q}_{N}$, with $x_{i}\in\\{0,1\\},\forall
i\in\\{1,\ldots,N\\}$, where the complementary vector $\bar{\bf x}$ is defined
such that $\forall x_{i}\in\\{0,1\\}$, $\bar{x}_{i}:=1-x_{i}$ , $\forall
i\in\\{1,\ldots,N\\}$. We assume that the sequences are circular, i.e.
periodic boundaries: $x_{N+i}=x_{i},\forall i\in\\{1,\ldots,N\\}$. From a
biological point of view in this binary representation each bit is a locus in
a gene located on the first strand. In this sense, our model mimics the
molecular evolution of circular viruses (Solé and Elena,, 2018; Tisza et al.,,
2020). For computational purposes, we model the molecular evolution by
selection dynamics as a two-step process. The first step involves the “access
to a new genotype” by mutation, whereas the second step determines which
mutated genotypes will be selected (or fixed). Let us remark that we do not
include recombination, because we are modelling asexual replication.

### Genetic mutation operations

Besides classical single locus mutations (bit-flip), chromosomal inversions
contribute to another level of genomic variation, which consists of an order
reversal in a segment of the genome. In this case the double-stranded
characteristic of DNA plays an important role because a molecular inversion is
not a mere permutation between sites of the sequence. Rather, it consists in a
permutation of the complementary strand, to then be exchanged with the main
strand within the segment where the inversion occurs (Griffiths et al.,, 2012,
p. 614). We model the chromosomal-like inversion as a two step procedure: a
conjugation operation, followed by a permutation operation.

We define the conjugation operation $\mathrm{\hat{C}}$ such that

$\mathrm{\hat{C}}:(x_{i},x_{i+1},\ldots,x_{j-1},x_{j})\longrightarrow(\bar{x}_{i},\bar{x}_{i+1},\ldots,\bar{x}_{j-1},\bar{x}_{j}),$

for $i,j\in\\{1,...,N\\}$. Note that, as the genome is circular, there is no
relation of order between $i$ and $j$ (i.e. $\mathrm{\hat{C}}$ is well-defined
even when $i>j$). The permutation operation $\mathrm{\hat{P}}$ is defined as

$\mathrm{\hat{P}}:(x_{i},x_{i+1},\ldots,x_{j-1},x_{j})\longrightarrow(x_{j},x_{j-1},\ldots,x_{i+1},x_{i}),$

for $i,j\in\\{1,...,N\\}$. Then, we have the (two-step) inversion operation:

$\mathrm{\hat{I}}\equiv\mathrm{\hat{C}}\circ\mathrm{\hat{P}}.$ (1)

Trivially, it can be verified that $\mathrm{\hat{C}}$ and $\mathrm{\hat{P}}$
commutes:
$\mathrm{\hat{C}}\circ\mathrm{\hat{P}}=\mathrm{\hat{P}}\circ\mathrm{\hat{C}}$.
Besides, Eq. 1 also defines a point mutation, i.e. when $i=j$, then
$\mathrm{\hat{I}}$ is a single bit-flip, equivalent to the biological
transversion-type point mutation (Griffiths et al.,, 2012, p. 555). These
operations can be easily translated into a model of computation capable of
simulating the evolutionary dynamics of a genetic-like Boolean sequence of
length $N$.

### NK fitness landscape

After the picturesque notion of “fitness landscape”, introduced in
evolutionary biology by Wright, (1932), adaptive evolution is often
conceptualised as hill-climbing walks on rugged landscapes. This metaphor even
became a “tool” to quantify (and think) many features about molecular
evolution (Gavrilets,, 2004; Obolski et al.,, 2018). A well tailored model to
simulate artificial molecular evolution is the NK fitness landscape
(Kauffman,, 1993, Ch. 2).

In the NK model, besides the genome length $N\in\mathbb{N}$, there is an
integer $K\in\mathbb{Z}_{(0,N-1)}$, which describes the epistatic interaction
between the sites in the genome and the contribution of each component to the
total fitness, which depends on its own value as well as the values of $K$
other sites. The fitness per bit is formally defined as:

$f_{i}:\\{0,1\\}^{K+1}\longrightarrow[0,1),\,\,\,\,\forall 1\leq i\leq N.$

Here $f_{i}(x_{i},x_{i_{1}},\ldots,x_{i_{K}})$ depends on the state of site
$x_{i}\in\\{0,1\\}$ and $K$ other sites $x_{i_{K}}\in\\{0,1\\}$. The $f_{i}$’s
are given by $N\cdot 2^{K+1}$ independent and identically distributed random
variables sampled from a given uniform probability distribution. Different
features of the NK model with more general distribution functions have been
analysed by, for example, Evans and Steinsaltz, (2002), Durrett and Limic,
(2003) and Limic and Pemantle, (2004).

The two frequently used neighborhood models are full disordered and
$K$-regular graphs: (i) The random neighborhood model, where $i$ and $K-1$
other sites are chosen at random according to a uniform distribution from
$\\{1,2,\ldots,N\\}$; (ii) The adjacent neighborhood model, where $i$ and $K$
other sites are successively ordered, i.e. $i,i+1,\ldots,i+K$ (each variable
modulo $N$ when using periodic boundary conditions). The total fitness
$f\in[0,1)$ for the sequence ${\bf x}\in\mathcal{Q}_{N}$ is then defined as:

$f({\bf
x}):=\frac{1}{N}\sum_{i=1}^{N}f_{i}(x_{i},x_{i_{1}},\ldots,x_{i_{K}}),$ (2)

where $\\{i_{1},\ldots,i_{K}\\}\subset\\{1,\ldots,i-1,i+1,\ldots,N\\}$.

The most important feature of the NK model is that the parameter $K$ “tunes”
the landscape ruggedness, that is the distribution of fitness local maximums,
ranging from non-epistatic interactions when $K=0$ (a Mount-Fuji-like
landscape with a single peak), to the full rugged (or random) landscape when
$K=N-1$. This allows us to analyse the time scales as a function of the
“disorder” of local maximums (fitness peaks) over the landscape.

### Adaptive walks

The evolutionary dynamics of a single digital organism can easily be
translated in terms of adaptive walks over the fitness landscape. That is,
once given a starting sequence ${\bf x}\in\mathcal{Q}_{N}$, it varies through
successive mutations (one bit-flip or one inversion per time step) resulting
in a mutated sequence ${\bf y}\in\mathcal{Q}_{N}$. Then, the fitness $f({\bf
x})$ and $f({\bf y})$ are calculated according to the NK model (2). If $f({\bf
y})>f({\bf x})$, the mutated sequence is selected, otherwise other mutations
on ${\bf x}$ are tested until the fitness increases. Therefore, fitness
decreasing mutations cannot be fixed, analogous to the “strong selection, weak
mutation” regime in population genetics (Gerrish and Lenski,, 1998). The
consequence of this constraint is that once a local maximum of fitness is
reached, no point mutation will allow the evolutionary dynamics to continue.
It is just stuck in a local peak. Unless a beneficial inversion-type mutation
is involved, and possibly opens a new way to increase the fitness. This would
be equivalent to “escape” from the local peak. In this new version of adaptive
walks on NK landscapes, we include the probability $p\in[0,1]$ associated with
the occurrence of inversion-like mutations between two sites $i$ and $j$ of
the genome, such that $i\neq j\in\\{1,\ldots,N\\}$ (and the probability of
point mutations is $q=1-p,\therefore i=j$). When $p=1$, all the mutations are
inversions (without point mutational events, i.e. $q=0$), and vice versa. If
$0<p<1/2$, the inversions are less frequent than the point mutational events.
Then, the mutation interval $[i,j]\in\\{1,\ldots,N\\}$ is drawn from a pseudo
random number generator. Each time-step is determined by a mutational (random)
event constrained to an incremental fitness selection. With these
computational recipes we can simulate successive mutated genotypes and their
fitness values, to generate fitness time series of adaptive evolutionary
“hikes” for different instances of the fitness landscapes and different
initial genomes.

#### Experimental setup

Figure 1: Hikes schematic. The vertical axis represents generations (i.e.
evolution time-steps). Each point indicates when a mutation is fixed. Green
dots represent “clones” of a given initial genome ${\bf x}_{0}$. Blue and red
dots symbolise point and inversion mutations respectively. Three time scales
are outlined. An initial transient period in which point mutations
predominate. Followed by a “stasis”-like period, where the occurrence of
mutations is very low or null. Finally, a period that “re-initiates” the
evolutionary dynamics, which is only possible if mutation by inversion occurs.
The idea of evolutionary burst is outlined, which is roughly a succession of
point mutations after an inversion occurs.

To elucidate the relationship between the evolutionary time scales and the two
types of mutations considered in our model, we implemented a series of
repetitions of the adaptive walk described above to create a pool of
trajectories over the fitness landscape. That is, from a given initial genome
${\bf x}_{0}\in\mathcal{Q}_{N}$ of length $N$ we run different instances of
the evolutionary process up to maximum number of time steps $\tau_{max}$.
Theses instances are like a branching random walk, but instead of several
descendants from a given genome in each generation (like in population
genetics), we generate $2N$-hikes starting at the same ${\bf x}_{0}$
(imitating independent evolutionary experiments with “clones” of a same
initial genome ${\bf x}_{0}$). Therefore, each hike is an independent
exploration of accessible paths over a (sub)domain of the fitness landscape.
These samples are quenched in the sense that we do not run simulations of
trajectories for a given a set $\\{{\bf x}_{0}\\}\in\mathcal{Q}_{N}$ of
independent initial genomes, but, as mentioned above, each time-step is driven
by a mutational (random) event. They are composed by two probabilistic
processes, the first one is a random choice of the type of mutation (bit-flip
or inversion) accordingly to the probability $p$; in the second one, the sites
to be mutated are raffled. For example, this can by done seeding a pseudo
random numbers generator independently for each process. See figure 1 for a
schematic of the time evolution of hikes.

## Results

Now we summarise the results of our empirical exploration. We display
illustrative results111Data and code are available at
https://gitlab.inria.fr/letrujil/NK_hikes. for genomes of length $N=50$.

##### Fitness benchmark

First, before exploring the time scales, we gauged the fitness values for
different instances of NK landscapes, in function of $K$ for adjacent and
random neighbourhoods. We analysed three mutational regimes: inversions and
point mutations only ($p=1$ and $p=0$ respectively), and the mixing of both
operations ($p=1/2$). In figure 2 we can verify that when $K=0$, the expected
fitness value is $\langle f\rangle\simeq 2/3$, and when $K$ increases, then
the expected fitness value decreases (for point mutations $\langle
f\rangle\rightarrow 0.60$ and for inversions $\langle f\rangle\rightarrow
0.63$). The trends displayed in figure 2 are congruent with the well known
“complexity catastrophe” in the NK model, i.e. the tendency of (averaged)
local fitness peaks to decrease when the epistatic interactions increase
(Kauffman,, 1993, p. 52). This trend depends on the probability distribution
for the contributions to fitness $f_{i}$ (which in our case is an uniform
distribution of identically distributed random variables) (Solow et al.,,
1999).

Figure 2: Fitness benchmark. Expected fitness value for genomes of length
$N=50$ versus the epistatic parameter $K=(0,\ldots,49)$. Blue dots correspond
to point mutations only ($p=0$) while the red ones correspond to inversions
only ($p=1$). Green dots are for equiprobable occurrence of inversions and
point mutations ($p=1/2$). Left: adjacent epistatic interactions. Right:
random epistatic interactions. Each point is the average for $100$ independent
hikes.

Figure 2 also shows that when the system evolves under the sole inversion
regime, fitness values are generally higher than in the point mutation regime.
However, this phenomenon is stronger in the random neighbourhood model: When
the epistatic interactions are fully random, the fitness values in the point
mutation regime are higher than those of the inversion regime below $K\simeq
25$. This is likely to be due to the fact that inversions are local operators,
hence being more efficient when epistatic interactions are local. In the
following of the paper all simulations and analyses use the adjacent epistatic
neighbourhood.

### Inversions trigger evolutionary bursts

##### The access to new fitness peaks

Following the experimental protocol described above, we evaluate $100$ hikes
initialized with a same random genome $x_{0}$ for $\tau_{max}=2\times 10^{5}$,
$N=50$ and $K=20$. Figure 3 (a) shows a benchmark of the evolutionary process
with point mutations alone. The total final time per hike (i.e. the time when
the last mutation has been fixed) doesn’t even reach $10^{3}$ generations,
verifying that the hikes stop at local fitness peaks (i.e. from $\sim 10^{2}$
generations up to $\tau_{max}$ nothing else happens). Moreover, the elapsed
times between mutations are generally of some tens of generations (see for
example the top left panel in fig. 4). When inversions are allowed (i.e.
$p\neq 0$, figure 3 (d)) we can verify that the evolutionary dynamics lasts
much longer than when $p=0$. This can be appreciated in figure 3 (b, e).
Inspecting this figure, we can notice that there are hikes that, like the
previous case, stay blocked in a local maximum of fitness. However, and it is
the key finding of the present work, there exist late mutational events
–driven by inversion mutations– allowing the evolutionary process to continue
far away on time. This can be better appreciated in figure 3 (b, e), where the
tails of the distribution of the final time per hike are larger when $p\neq 0$
than when $p=0$. The total number of generations when inversions are fixed can
be almost $30$ times longer than the total number of generations reached when
point mutations are the only source of genome variation in the evolutionary
dynamics. Moreover, the action of inversions can drive the dynamics to higher
fitness values. This is shown in figure 3 (f), where the maximum value of
fitness in the case $p\neq 0$ is higher than the case $p=0$ (figure 3 (c)).
This shows that some “new” fitness peaks are “discovered” in comparison with
the hikes when $p=0$. In the histogram this corresponds to the bins of fitness
between $0.68$ and $0.73$. So we can paraphrase that not only can evolution
“go further”, it can also discover something new: Peaks that were previously
inaccessible over the fitness landscape, are now accessible thanks to the
action of inversions that provoke large-scale reorganisations of the genome of
the digital organism simulated in this work.

Figure 3: Effect of inversions on the evolutionary dynamics for $N=50$,
$K=20$, $\tau_{max}=2\times 10^{5}$ generations for a given initial genome
${\bf x}_{0}$. Top: $p=0$ (bit-flip only); Bottom: $p=0.1$ ($10$% of the
mutational events are inversions). (a) and (d) are scatter plots of the
mutational events per hike and for the two mutational conditions (blue open
square: bit-flips; red open circles: inversions). The histograms (b) and (e)
count the distribution of the final time step (i.e. generation of the last
fixed mutational event, being it a bit-flip or an inversion) for the two
mutational conditions. (c) and (f) display the distribution of the final
fitness values at generation $2\times 10^{5}$ for the two mutational
conditions. The red lines indicate the respective mean values.  Figure 4:
Fitness bursts. Fitness time-series of two representative hikes for $N=50$,
$K=20$ and $\tau_{max}=2\times 10^{5}$. The pale blue rows emphasise the
genomes muted by a bit-flip (highlighted with pale pink), while the pale pink
rows emphasise inversion mutational events. The underlined mutated segments in
the genome are framed by brackets, and the bits are boldfaced. Top left panel:
shows the generations when the mutational events were fixed for hike 65, when
$p=0.0$. Right panel: illustrates the case when $p=0.1$ for hike 54, showing
the stasis-like period and the occurrence of the evolutionary bursts. Bottom
left figure: Displays the fitness time-series of hike 54. See main text, and
figure 1, for the description of the $4$ evolutionary phases.

##### An explicit example of stasis and bursts

In figure 4 (bottom left and right panels) we display an example from the pool
of hikes when $p=0.1$ (for hike $54$ in figure 3 (d)). It can be seen that the
exploration of the fitness landscape is initially carried out by point
mutations until the sixty-sixth generation (this time-domain is labeled as $1$
in the figure). Then the dynamics remain in an invariant state for a
relatively long period of time that lasts until the time-step $789$ (time-
domain $2$). At that time an inversion is fixed which is rapidly followed by
the fixation of a point mutation at generation $899$. Then, the evolutionary
process returns to a “stasis-state” for a long period (between generations
$899$ and $28,735$ – time-domain $3$). Finally, at generation $28,735$ an
inversion occurs that boosts evolution again (time-domain $4$). A very
interesting observation is that at this stage there are a series of successive
point mutations (similar to the number of point mutations that occurred at the
beginning of evolutionary dynamics). We call this sudden mutational events an
evolutionary burst. This can be interpreted, using the fitness landscape
metaphor, as if the system has “jumped” from a local fitness peak to a new
one, hence discovering a new domain of the landscape. Then, the system begins
to explore and climb this new peak until it reaches a local (with respect to
Hamming distance) optimum of fitness. In this example, once again a new
inversion occurs (at generation $36,073$) and the fitness keeps increasing
again for a while. As we can see in figure 4, bottom-left panel, if we compare
with the early stages of the evolutionary process (time-domain $1$), the gain
of fitness is larger on domain $4$ than on the first one where only point
mutations were fixed. This shows that, when an inversion take place, it can
boost evolution.

### Effect of the epistasis degree K

Figure 5: Effect of the degree of epistasis. Scatter plots of the mutational
events per hike for $N=50$, $p=0.01$, $\tau_{max}=2\times 10^{6}$ generations.
The simulated instances, from “low” to “high” landscape-ruggedness, are: (a)
$K=2$, (b) $K=10$, (c) $K=20$, (d) $K=30$, (e) $K=40$, (f) $K=49$.

##### From low to high ruggedness

In figure 5 we compare the behaviour of the evolutionary dynamics for the NK
model with different values of the epistatic parameter
$K=\\{2,10,20,30,40,49\\}$ ($N$ still being fixed to $50$). That is, from a
moderately rugged landscape to a full random one. In this case the total
number of hikes is $100$ (all starting with a simple genome made of a string
of $50$ zeros), $p=0.01$ and $\tau_{max}=2\times 10^{6}$ generations. In this
analysis, inversions occur once every hundred times in average, allowing us to
better appreciate the effects of the fitness landscape ruggedness on
timescales and the duration of the state of stasis due to the epistatic
interaction between the $K$ loci. For the epistatic interaction values
analysed here, it can be seen that during the first stages of the evolutionary
process, point mutations always drive evolution, corresponding to an efficient
initial exploration of accessible paths (by “efficient” we mean that point
mutations can be easily and rapidly explored as there are only $N=50$
different possibilities). This period lasts roughly from $10^{2}$ to $10^{3}$
generations, even though this is less differentiated for the full random
landscape ($K=49$). Then, for $K=2$, $10$ and $20$, we can identify a gap of
evolutionary activity (the stasis-like period) that lasts approximately up to
$10^{4}$ generations. Some mutations occur during this period, but most of
them were preceded by inversions that took place early in the evolutionary
process. Then, from $K=$ 20, the mutations are more scattered and the stasis-
like period tends to vanish. When $K$ is large, the number of local fitness
peaks increases; therefore, even when genomes are able to escape their current
peak, it is likely they have only found a different local optima rather than a
global optima. Therefore, since inversions only release genomes from local
optima rather than finding global optima, such inversion mutations are less
evolutionarily efficient as $K$ increases. $K$ also modulates the number of
mutational events that occurs during an evolutionary burst (including the
initial one), less mutational events being fixed when $K$ increases.

##### Examples of time scale dependence with epistasis

Finally, in figure 6 we present a random selection, from the experiments
presented above, of those events that can be identified as possible
evolutionary bursts. Here we can see more clearly how inversion mutational
events contribute to escaping local fitness peaks, especially in the case when
$K=2$, where it is noted that the different hikes quickly get stuck in local
fitness peaks: as the landscape is smoother point mutational events are fixed
more often. Once again we see how the evolutionary dynamics reaches a state of
stasis, evolution being boosted again with the occurrence of an inversion
after a long period of waiting time. For $K=10$, $20$ and $30$, the magnitudes
of the fitness increases are quite large, in comparison to the case when
$K=2$, which translates into a highly steep slope, giving the impression that
the evolutionary process has not only been re-started but also boosted. For
the values of $K=40$ and $49$, the stasis period has almost completely
disappeared, and mutational inversion events occur at any evolutionary time
scale. This is easily understood, since the increase in the disorder of the
(now many) fitness peaks has an effect on the “randomization” of evolutionary
steps. When epistasis is high, evolutionary dynamics is more “turbulent”,
meaning that there are no well-differentiated time scales between the types of
mutations. We can also see in the scatter plots of the mutational events
(figure 5) the progressive “melt” of differentiated time scales, where the
duration of the stasis is also appreciated as a function of the degree of
epistatic interactions.

Figure 6: Stasis and burst. Fitness time-series for a random selection of
hikes picked from the simulations reported in fig. 5, showing the stasis- and
burst-like events: (a) $K=2$, (b) $K=10$, (c) $K=20$, (d) $K=30$, (e) $K=40$,
(f) $K=49$.

## Conclusions

From the empirical analysis of our model, it has been shown that complex long-
and short-term evolutionary dynamics can be simulated on rugged fitness
landscapes. To achieve that, we have extended the well-known family of NK
landscapes by including inversion-like chromosome rearrangements (in this
sense, our work complements somewhat related methods to simulate the NK model
beyond bit-flips mutations alone, e.g. (Hadany and Beker,, 2003; Aguirre and
Tanaka,, 2003; Østman et al.,, 2012; Gillespie et al.,, 2018)). We then
studied the effects of epistasis on the evolutionary time scales of digital
genomes.

The inversion-like mutations have been incorporated into the simulation of the
evolutionary dynamics of NK digital genomes through the probability $p$, which
controls the relative frequency of occurrence of each kind of mutation. By
combining the effect of $p$ and $K$ (level of epistasis), we have been able to
reproduce behaviors similar to long waiting periods in evolutionary dynamics –
suggestively evoking evolutionary stasis – followed by mutational bursts.
Therefore, with this in hand we have been able to model the length of adaptive
paths as a composition between two different mutational events, the
traditional bit-flip that explores the neighborhood of a short-range domain of
the fitness landscape, and the effects of genome sequence inversion, akin to
long-range exploration (Bagnoli and Bezzi,, 2001; Tomassini,, 2016). All this,
together with the combinatorial aspect of all possible variations of any
binary sequence, also determine the orders of magnitude of the time scales.
Indeed, combinatorics imposes its time scale which depends on the size of the
binary sequence and on the specifics of the different mutational operations
and therefore deserves to be studied in more detail than the results presented
here.

At least three main time scales were found in our simulations. The first one,
is an exploration mainly driven by bit-flips up to a local maximum of fitness.
Then, a stasis-like regime where almost no mutation is fixed. This stationary
state can be broken with the fixation of an inversion-like mutation, allowing
for the exploration of a new region of the fitness landscape. Once the digital
organism “jumps off” a local peak, it can reach a new “hill” in the fitness
landscape. This opens the possibility of a new exploration, which is again
carried out mainly by point mutations that may lead to burst-like events, i.e.
the rapid succession of point mutations following an inversion. Indeed, in the
present model bursts can only be triggered by chromosomal inversions but real
biological sequences may undergo a vast repertoire of chromosomal
rearrangements (inversions, translocations, duplications…) that deserves
future studies. The system may then alternates between phases of stasis and
bursts (triggered by inversion like mutations) up to a “dead-end” where all
accessible favorable paths have been followed and the evolutionary potential
of the (last mutated and selected) binary genome has been exhausted and/or
blocked by epistatic interactions. Note that such “dead-end” may be relatively
quickly reached in our simulations (the size of the mutational repertoire
being proportional on the genome size $N$) but that, in longer/real genomes,
the number of possible inversions ($N^{2}$) quickly grows. We suggest that
this enables virtually infinite exploration of new possibilities, akin to
open-endedness (Banzhaf et al.,, 2016). This is another topic that deserves to
be investigated and that is also a subject of debate in the evolutionary
biology and artificial life communities.

The fitness time-series reported in this work are examples of the statistics
of records increments, and the coexistence of short- and long-range time
scales resembles the drift and diffusion jumps phenomenology. That is a
Brownian-like random walk including Lévy flights (Masuda et al.,, 2017). This
was also evidenced in the scatter plots of figure 5 and the transition-like
behaviour towards a more “turbulent” evolutionary dynamic modulated through
epistasis.

Finally, our model reproduces an evolutionary pattern alternating short and
long time scales, without resorting to a description at the population level.
It is often supposed that the discrepancy between short and long evolutionary
time scales is due to long-term environmental variations (Uyeda et al.,,
2011). Despite its simplicity, our model shows that such patterns can occur in
stable environments through the sole intricacy of different mutation
operations having different exploration time scales. The present work
highlight that mutations which locally reverse the genome order mimic the
process called by Kauffman, (2000) as “expanding the adjacent possible”. In
our case, this can be translated as “inversion-like mutation driving digital
organisms into a new domain on the fitness landscape”, hence ending possibly
long-lasting stasis periods and triggering bursts of bit-flip mutations. In
this sense, our model could describe a simple scenario that, nonetheless, can
shed new light on punctuated equilibrium and complex interlocking micro-
macroevolutionary dynamics.

## References

  * Aguirre and Tanaka, (2003) Aguirre, H. E. and Tanaka, K. (2003).  Genetic algorithms on nk-landscapes: Effects of selection, drift, mutation, and recombination.  In Workshops on Applications of Evolutionary Computation, pages 131–142. Springer.
  * Bagnoli and Bezzi, (2001) Bagnoli, F. and Bezzi, M. (2001).  Small world effects in evolution.  Physical Review E, 64(2):021914.
  * Banzhaf et al., (2016) Banzhaf, W., Baumgaertner, B., Beslon, G., Doursat, R., Foster, J. A., McMullin, B., De Melo, V. V., Miconi, T., Spector, L., Stepney, S., and White, R. (2016).  Defining and simulating open-ended novelty: requirements, guidelines, and challenges.  Theory in Biosciences, 135(3):131–161.
  * Casanova and Konkel, (2020) Casanova, E. L. and Konkel, M. K. (2020).  The developmental gene hypothesis for punctuated equilibrium: combined roles of developmental regulatory genes and transposable elements.  Bioessays, 42(2):1900173.
  * Darwin, (1964) Darwin, C. (1964).  On the origin of species: A facsimile of the first edition, volume 11.  Harvard University Press.
  * Durrett and Limic, (2003) Durrett, R. and Limic, V. (2003).  Rigorous results for the nk model.  Annals of probability, 31(4):1713–1753.
  * Evans and Steinsaltz, (2002) Evans, S. N. and Steinsaltz, D. (2002).  Estimating some features of $nk$ fitness landscapes.  Annals of Applied Probability, 12(4):1299–1321.
  * Gavrilets, (2004) Gavrilets, S. (2004).  Fitness landscapes and the origin of species (MPB-41).  Princeton University Press.
  * Gerrish and Lenski, (1998) Gerrish, P. J. and Lenski, R. E. (1998).  The fate of competing beneficial mutations in an asexual population.  Genetica, 102:127–144.
  * Gillespie et al., (2018) Gillespie, L., Dolson, E., Lalejini, A., and Ofria, C. (2018).  Changing environments drive the separation of genes and increased evolvability in nk-inspired landscapes.  In The 2018 Conference on Artificial Life.
  * Gould and Eldredge, (1977) Gould, S. J. and Eldredge, N. (1977).  Punctuated equilibria: the tempo and mode of evolution reconsidered.  Paleobiology, pages 115–151.
  * Griffiths et al., (2012) Griffiths, A. J. F., Wessler, S. R., Carroll, S. B., and Doebley, J. (2012).  Introduction to genetic analysis.  W. H. Freeman.
  * Hadany and Beker, (2003) Hadany, L. and Beker, T. (2003).  Fitness-associated recombination on rugged adaptive landscapes.  Journal of evolutionary biology, 16(5):862–870.
  * Ispolatov et al., (2019) Ispolatov, I., Alekseeva, E., and Doebeli, M. (2019).  Competition-driven evolution of organismal complexity.  PLoS computational biology, 15(10):e1007388.
  * Kauffman and Levin, (1987) Kauffman, S. and Levin, S. (1987).  Towards a general theory of adaptive walks on rugged landscapes.  Journal of Theoretical Biology, 128(1):11–45.
  * Kauffman, (1993) Kauffman, S. A. (1993).  The origins of order: Self-organization and selection in evolution.  Oxford University Press, USA.
  * Kauffman, (2000) Kauffman, S. A. (2000).  Investigations.  Oxford University Press.
  * Kauffman and Weinberger, (1989) Kauffman, S. A. and Weinberger, E. D. (1989).  The nk model of rugged fitness landscapes and its application to maturation of the immune response.  Journal of theoretical biology, 141(2):211–245.
  * Limic and Pemantle, (2004) Limic, V. and Pemantle, R. (2004).  More rigorous results on the kauffman–levin model of evolution.  The Annals of Probability, 32(3):2149–2178.
  * Masuda et al., (2017) Masuda, N., Porter, M. A., and Lambiotte, R. (2017).  Random walks and diffusion on networks.  Physics reports, 716:1–58.
  * Obolski et al., (2018) Obolski, U., Ram, Y., and Hadany, L. (2018).  Key issues review: evolution on rugged adaptive landscapes.  Reports on Progress in Physics, 81(1):012602.
  * Østman et al., (2012) Østman, B., Hintze, A., and Adami, C. (2012).  Impact of epistasis and pleiotropy on evolutionary adaptation.  Proceedings of the Royal Society B: Biological Sciences, 279(1727):247–256.
  * Pennell et al., (2014) Pennell, M. W., Harmon, L. J., and Uyeda, J. C. (2014).  Is there room for punctuated equilibrium in macroevolution?  Trends in ecology & evolution, 29(1):23–32.
  * Solé and Elena, (2018) Solé, R. and Elena, S. F. (2018).  Viruses as complex adaptive systems, volume 15.  Princeton University Press.
  * Solow et al., (1999) Solow, D., Burnetas, A., Tsai, M.-C., and Greenspan, N. S. (1999).  Understanding and attenuating the complexity catastrophe in kauffman’s nk model of genome evolution.  Complexity, 5(1):53–66.
  * Tisza et al., (2020) Tisza, M. J., Pastrana, D. V., Welch, N. L., Stewart, B., Peretti, A., Starrett, G. J., Pang, Y.-Y. S., Krishnamurthy, S. R., Pesavento, P. A., McDermott, D. H., et al. (2020).  Discovery of several thousand highly diverse circular dna viruses.  Elife, 9:e51971.
  * Tomassini, (2016) Tomassini, M. (2016).  Lévy flights in neutral fitness landscapes.  Physica A: Statistical Mechanics and its Applications, 448:163–171.
  * Uyeda et al., (2011) Uyeda, J. C., Hansen, T. F., Arnold, S. J., and Pienaar, J. (2011).  The million-year wait for macroevolutionary bursts.  PNAS, 108(38):15908–15913.
  * Wright, (1932) Wright, S. (1932).  The roles of mutation, inbreeding, crossbreeding, and selection in evolution.  In Proc. Sixth International Congress on Genetics, volume 1, pages 356–366. na.


# Algorithmic Recourse for Dependent Data

Anonymous Authors  Email AffiliationsCountries

(2018)

###### Abstract.

Algorithmic Recourse, Counterfactual Explanations, Explainable Machine
Learning

††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock,
NY††booktitle: Woodstock ’18: ACM Symposium on Neural Gaze Detection, June
03–05, 2018, Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs:
Computing methodologies Learning latent representations††ccs: Computing
methodologies Knowledge representation and reasoning

##  1\. INTRODUCTION

##  2\. RELATED WORK

##  3\. BACKGROUND

##  4\. Our Framework: DEAR

##  5\. CONCLUSION