#  gwfish: A simulation software to evaluate parameter-estimation capabilities
of gravitational-wave detector networks

Jan Harms1,2 Ulyana Dupletsa1,2 Biswajit Banerjee1,2 Marica Branchesi1,2 Boris
Goncharov1,2 Andrea Maselli1,2 Ana Carolina Silva Oliveira3 Samuele
Ronchini1,2 Jacopo Tissino1,2 1Gran Sasso Science Institute (GSSI), I-67100
L’Aquila, Italy  2INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi,
Italy  3Department of Physics, Columbia University in the City of New York,
New York, NY 10027, USA

###### Abstract

An important step in the planning of future gravitational-wave (GW) detectors
and of the networks they will form is the estimation of their detection and
parameter-estimation capabilities, which is the basis of science-case studies.
Several future GW detectors have been proposed or are under development, which
might also operate and observe in parallel. These detectors include
terrestrial, lunar, and space-borne detectors. In this paper, we present
gwfish111github.com/janosch314/GWFish, a new software to simulate GW detector
networks and to calculate measurement uncertainties based on the Fisher-matrix
approximation. gwfish models the impact of detector motion on PE and makes it
possible to analyze multiband scenarios, i.e., observation of a GW signal by
different detectors in different frequency bands. We showcase a few examples
for the Einstein Telescope (ET) including the sky-localization of binary
neutron stars, and ET’s capability to measure the polarization of GWs.

##  I Introduction

The gravitational-wave (GW) community is currently in the phase of developing
science cases and analysis tools of future GW detector networks [1, 2]
including potential upgrades of the current infrastructures Virgo [3], LIGO
[4], LIGO India [5], KAGRA [6], and the proposed Einstein Telescope [7, 8] and
Cosmic Explorer [9]. The approved space-borne detector LISA is expected to
begin observations in the second half of the 2030s [10]. Entirely new detector
concepts are under study like the Lunar GW Antenna (LGWA) [11]. The networks
can be formed between alike detectors observing at the same frequencies, but
also combining observations in different frequency bands of the same GW
sources [12, 13, 14]. Providing a simulation of conceivable observation
scenarios formed by these detectors is an important part of the science-case
development of future GW detector networks, but also a challenging and
potentially computationally expensive task. In recent years, a few so-called
Fisher-matrix codes were developed for this purpose [15, 14, 16].

In this article, we present the simulation software gwfish, which uses
frequency-domain GW models combined with a time-domain simulation of GW
detector networks through the stationary-phase approximation. This creates a
framework to study important aspects of PE with future detector networks,
where the change of position and orientation of a detector during the
observation of a signal can have an important impact especially on the sky
localization [17, 15, 14, 18]. The calculation of PE errors is done using the
Fisher-matrix approximation [19], which corresponds to a Gaussian
approximation of the likelihood expected to be acceptable for signals with
higher signal-to-noise ratio. When the Fisher matrix is used to directly
evaluate PE errors like in gwfish, then priors cannot be considered. As a
consequence, overestimation and underestimation of PE errors are possible
[20]. It should be noted that Gaussian likelihood approximations based on
Fisher matrices can also be used for posterior sampling, which leads to a
significant speed-up of the likelihood evaluation. This method can be combined
with arbitrary priors.

A feature of gwfish is that multiband scenarios can be simulated as well. It
is a simple step to do with Fisher-matrix codes since Fisher matrices for
different detectors are added irrespective of the frequency band that provided
the signal information. Technically, it just requires to set up the code so
that detector motion and signal waveforms can be simulated accurately in
different frequency bands. This includes the application of the time-delay
interferometry (TDI) formalism to simulate space-borne detectors like LISA
beyond the long-wavelength regime [21, 22, 23], and the support of multiple
components of a detector, e.g., multiple interferometers of a xylophone
configuration [24] or multiple sensors for LGWA type detectors [11]. Each
component can be assigned a duty cycle for more realistic assessments of the
observing scenario. Details of the detector simulation are described in
section II.

The Fisher-matrix formalism is computationally efficient, but numerically less
robust than posterior sampling. The numerical challenges concern the
calculation of waveform derivatives and the inversion of Fisher matrices.
Waveform derivatives are best carried out in a hybrid analytical-numerical
scheme to minimize the computational effort and numerical errors, and to keep
it compatible with arbitrary waveform models (also referred to as waveform
approximants). In gwfish, the numerical differentiation is tuned to the
properties of waveform models to reduce numerical errors. Concerning the
inversion of Fisher matrices, the main issue is that these matrices can be
very close to singular, i.e., with a huge range of possible eigenvalues, which
is connected to signal-model degeneracies. To get accurate PE errors at least
of the parameters not involved in the model degeneracy, one needs to
systematically deal with matrix singularity. The general PE framework of
gwfish and how it addresses the numerical challenges are described in section
III.

As a demonstration, we show basic PE results for Einstein Telescope (ET) in
section IV, which is a good example to understand the capabilities of gwfish,
but also its limitations since single-detector scenarios are especially
difficult to model in the Fisher-matrix formalism due to unavoidable
degeneracies in the signal model. As a direct scientific conclusion of these
analyses, we confirm that the ET configuration has unique capabilities for PE
of compact-binary coalescences. Other more detailed studies carried out with
gwfish have been published in separate papers. These include the evaluation of
the perspectives of multi-messenger observations for ET (as single observatory
and in network of GW detectors) operating in synergy with $\gamma$ and X-ray
satellites using astrophysically motivated populations of BNS mergers [25],
and very high-energy observatories (Banerjee et al. 2022 in preparation).
gwfish was also used for the simulation of an astrophysical background in the
ET null-stream study [26].

##  II Detector-network simulation

The simulation of a detector network needs to consider several aspects, which
are represented under different categories inside gwfish:

  * $\bullet$

Component (laser interferometer, TDI channel, accelerometer): time-varying
detector response to GWs, noise modeling, duty cycle;

  * $\bullet$

Detector: one or multiple components observing the same frequency band and
added as a group to a detector network;

  * $\bullet$

Network: one or multiple detectors observing together.

Most of the effort of the detector simulation goes into the calculation of the
detector response. We outline the basics and refer to other publications for
details. The long-wavelength approximation is applied to all terrestrial GW
detectors simulated with gwfish, which assumes that the arm-length of a GW
detector is much smaller than the length of a GW. It is not a fully accurate
approximation above a kilohertz for a 40 km detector like the proposed Cosmic
Explorer [27, 28, 29]. The response formalism for space-borne detectors can be
adapted to ground-based detectors, e.g., to calculate the GW response of a
Fabry-Perot cavity beyond the long-wavelength regime [27].

The response of a detector is complicated by the fact that the detector
orientation and position can change with time. Especially the change in
orientation is important to simulate. Since signals are represented in
frequency domain in gwfish, the question is how to convolve the detector
response with the signal. In gwfish, it is assumed that a function $t(f)$
exists mapping signal frequencies to times. To obtain this function, we
exploit the stationary phase approximation, which underlies most analytical
models of compact-binary signals and is known to lead to very accurate
approximations of the Fourier spectrum of compact-binary signals [30]. This
function can be calculated from the phase $\phi(f)$ of the waveform:

$t(f)=\frac{1}{2\pi}\frac{{\rm d}\phi(f)}{{\rm d}f}.$ (1)

A technical issue with this approach to calculate $t(f)$ is that it requires a
dense sampling in frequency space to have enough resolution of the phase
$\phi(f)$ for a finite-difference method. Since we do not want to impose any
unnecessary constraints on waveform modeling in gwfish, we instead calculate
$t(f)$ directly from the lowest-order term of the phase evolution, e.g., see
equation (15) in [31]. This approximation is acceptable since only the early
phase of the inspiral needs to be modeled accurately to account for the
effects of detector motion on PE. Once $t(f)$ is calculated for a signal, it
can be used as a time vector for the simulation of the motion of the
detectors. The GW signal of a detector can then be written in terms of a
response tensor $\mathcal{A}(t(f))$ contracted with the GW strain tensor
$h(f)$ [32]

$\mathcal{R}(f)=\mathcal{A}(t(f)):h(f).$ (2)

For example, the response tensor for a laser interferometer with two arms
along direction $\vec{e}_{1},\,\vec{e}_{2}$ assumes the form

$\mathcal{A}=\frac{1}{2}(\vec{e}_{1}\otimes\vec{e}_{1}-\vec{e}_{2}\otimes\vec{e}_{2}),$
(3)

where all quantities are functions of $t(f)$. It should be noted that in the
case of terrestrial detectors, the actual frequency-dependent response of a GW
detector, e.g., any aspect going into detector calibration, is included in the
instrument-noise spectrum referred to GW strain.

Figure 1: Detectors included in gwfish with simulated detector motions
indicated by arrows.

For LISA-type space detectors, another method must be used to calculate the GW
response since the long-wavelength approximation does not hold [21]. The
fundamental building blocks are the readouts $y_{ij}(f)$ of individual
spacecraft links $j\rightarrow i$, which are then electronically combined with
time delays to mimic the workings of a laser interferometer. The main purpose
of combining different readouts is to reduce laser-frequency noise by several
orders of magnitude, which is achieved in terrestrial detectors by optical
interference. This method is called time-delay interferometry (TDI).
Effectively, TDI makes each satellite of a triangular detector configuration a
vertex of a laser interferometer. It turns out that the instrument noise in
the basic Michelson-type TDI combinations is correlated, which means that it
is convenient to diagonalize the noise-correlation matrix, which yields the
three (noise) independent TDI channels $A,\,E,\,T$ named after the three
scientists who set up the TDI formalism for LISA (Armstrong, Estabrook, and
Tinto). The $T$ channel has a strongly suppressed response to GWs below
$f_{0}=c/(2L)$, where $L$ is the average arm length of the space detector. The
channel $T$ is the sum of the three Michelson-type TDI channels, and it is
known as null stream in terrestrial GW detectors [33, 34]. The three
independent channels $A,\,E,\,T$ are simulated in gwfish for space detectors.
An approximation made in the calculation is that the so-called breathing
motion of the space-detector configuration, i.e., a small modulation of the
arm lengths over the course of a year, is neglected. Correcting the TDI
combination for this effect is important for the suppression of laser-
frequency noise, but it has a small effect on signal detection and PE. In
other words, the satellite configuration is simulated as rigid body orbiting
the Sun and with a rotating detector plane known as cartwheel motion [23].

The third detector type simulated in gwfish is LGWA, which consists of an
array of seismometers measuring the surface displacement of the Moon produced
by a passing GW. Here, the long-wavelength approximation applies again since
the detector is designed for measurements in the decihertz band and below, and
the length of the GWs need to be compared with the diameter of the Moon.
Therefore, the response of a seismometer to GWs can again be described by
equation (2), and the response tensor takes the form

$\mathcal{A}=\vec{e}_{n}\otimes\vec{e}_{1},$ (4)

where $\vec{e}_{1}$ is the displacement direction measured by the seismometer,
and $\vec{e}_{n}$ is the surface normal vector at the seismometer position.
The time dependence of the response comes from the rotation of the Moon. The
effect of its orbital motion around Earth is neglected.

##  III Parameter estimation

The PE is simulated in gwfish through the Fisher-matrix formalism, which is
based on a Gaussian approximation of the likelihood function:

$\mathcal{L}(\vec{x})\propto\exp(-\vec{x}^{\,\intercal}\cdot\mathcal{C}^{-1}\cdot\vec{x}/2).$
(5)

The covariance matrix $\mathcal{C}$ is the inverse of the Fisher matrix
$\mathcal{F}$, and its components can be calculated from the signal model
$h(\vec{x}\,)$,

$\mathcal{F}_{ij}=\sum\limits_{k=1}^{N}\langle\partial_{x_{i}}h^{k}|\partial_{x_{j}}h^{k}\rangle,$
(6)

where $N$ is the number of components of the GW detector network, and the
derivatives of the waveform are calculated with respect to the model
parameters $\vec{x}=(x_{i})$. The expression makes use of an inner product
defined on the signal-model manifold as follows:

$\langle a|b\rangle=4\int\limits_{0}^{\infty}{\rm
d}f\frac{\Re(a(\vec{x},f)b^{*}(\vec{x},f))}{S_{\text{n}}(f)}.$ (7)

The product requires a model of the instrument noise spectral density
$S_{\text{n}}(f)$. The theoretical framework is simple, but there are
practical challenges. First, for some of the model parameters, the derivatives
need to be calculated numerically. It is not necessary for all waveform
models, but one ought to adopt an approach that is valid for all models.
Numerical differentiation suffers from limitations of numerical precision,
e.g., some derivatives might be calculated from tiny variations of potentially
very large numbers. In gwfish, derivatives with respect to waveform phase at
merger $\phi_{\text{c}}$, merger time $t_{\text{c}}$, and source luminosity
distance $D$ are calculated analytically. Derivatives of all other parameters
(masses, spins, tidal deformation, sky location, inclination angle,
polarization angle) are calculated numerically. The waveform phase is a sum of
terms with large and small values. Especially the component $2\pi
ft_{\text{c}}$ can take very large numerical values, and gwfish eliminates
this term in calculations of numerical derivatives involving the phase of the
waveform.

Another numerical challenge is connected to the inversion of the Fisher matrix
in equation (6) to obtain the PE errors. The first issue is that model
parameters can have a large range of numerical values, e.g., $10^{30}$ for
masses in kg and values of order 1 for all angular variables. The situation
can be ameliorated by normalizing the Fisher matrix [35]:

  1. 1.

Collect the diagonal elements $d_{i}=\mathcal{F}_{ii}$ of the Fisher matrix.

  2. 2.

Divide row $k$ of the Fisher matrix by $\sqrt{d_{k}}$. Divide column $k$ of
the Fisher matrix by $\sqrt{d_{k}}$. This leaves a matrix with diagonal
components equal to one and all off-diagonal elements in the range $[-1,1]$.

  3. 3.

Invert the matrix, then repeat step (2) with the inverted matrix using the
same values for $d_{k}$.

A remaining potential numerical issue is when Fisher matrices are close to
degenerate. To address this issue, gwfish excludes all singular values lying
below a certain numeric precision threshold from the matrix before inverting
it [36]. The main outcome of this step is to make the matrix inversion robust
with respect to all significant singular values containing physical
information. In cases where components of small singular values need to be
discarded, the waveform model is close to degenerate with respect to at least
one parameter pair, e.g., distance and inclination angle can show strong
correlations. The PE errors estimated for such strongly correlated parameters
are typically substantially different from PE errors coming out of a full
Bayesian analysis (posterior sampling). For example, usage of prior
distributions in Bayesian analyses can substantially improve parameter
estimates in these cases. Model degeneracy is one of the fundamental failure
modes of the Fisher-matrix analysis.

gwfish can use all waveform models of LALSimulation [37], and it provides
independent implementations of the TaylorF2 (of 3.5th post-Newtonian order)
and IMRPhenomD approximants. It is beyond the scope of this paper to provide a
useful guideline for approximant usage with gwfish. Generally, it is important
to understand the physics that the LALSimulation approximant incorporate, and
these can have an enormous impact on parameter estimation. For example, it
should be clear that for the analysis of high-mass or high-redshift BBH
signals, which appear only with several cycles in the observation band, an
approximant like TaylorF2 is inappropriate since it does not include the
merger and ringdown of the binary. In these cases, the safest choice is to use
approximants that include higher-order modes. However, there are good use
cases even for TaylorF2, e.g., if one wants to calculate signal-to-noise ratio
(SNR) and localization of binary neutron stars. It is recommended to make
comparison of results for a substantial random sample of signals (of order 10
– 100 signals) before deciding, whether it is safe to use a simplified
waveform model.

###  III.1 bilby versus gwfish

We compare SNRs and PE between bilby — a detector simulation and Bayesian
analysis software [38] — and gwfish for four signals observed with ET (see
Table 1). Due to current constraints imposed by bilby, and by other state-of-
the-art PE software, it is not possible to carry out accurate Bayesian
analysis of BNS signals for ET. This has to do with computational limitations.
Standard representations of BNS waveforms for ET would require about 108
samples due to the day-long observation time of these signals with ET, and
bilby also does not consider Earth’s rotation during the observation of GW
signals for waveform modeling (an effect that does not play a role for
analyses of compact binaries with current GW detectors). As a consequence, we
carry out our comparison for BBHs, which, in some sense, is also more
interesting since model degeneracies play a more important role. We do the
comparison for each signal with a random realization of the detector noise in
bilby in order to compare the gwfish results with a more realistic simulation
of the Bayesian analysis. However, while noise realizations can generally
influence PE errors [39], it turns out that it has a negligible effect here
since the SNRs of all signals are very high.

Signals 1 – 3 have randomly sampled angular parameters, while the fourth
signal (Signal 3*) has the same parameters as Signal 3 except for the two
masses which are the same as for Signal 1. The signals are at relatively low
redshift to make sure that the SNRs are high. In this case, the Gaussian
approximation of the likelihood function is expected to hold, but we will see,
not entirely surprisingly, that PE errors produced by gwfish and bilby are
still not identical. As for the detector, we employ ET alone with its three
components. For both bilby and gwfish we are using the inspiral-merger-
ringdown waveform IMRPhenomD from LALSimulation, setting the minimum frequency
to 10 Hz and the maximum frequency to 1024 Hz. For gwfish, we include results
with 2 Hz lower bound to show the impact of ET’s low-frequency sensitivity on
PE. The reason not to use a 2 Hz lower bound for waveforms in bilby is that
bilby enforces a frequency resolution of about 1/(duration of signal), which
means that the number of samples used for each waveform becomes very large
even for BBHs with a 2 Hz lower bound compared to 10 Hz. Consequently, it
takes too long to finish posterior sampling. In our analysis, we used a
frequency resolution of 1/153 Hz for Signal 1, 2 and 3, and 1/108 Hz for
Signal 3*. A simpler waveform model such as TaylorF2 would have been
inappropriate here since model degeneracies would be higher especially when
simulating BBH observations involving only ET.

Table 1: Parameters of the signals used for SNR and PE comparisons between bilby and gwfish. The spins are assumed to be aligned with the orbital angular momentum. | Signal 1 | Signal 2 | Signal 3, 3*
---|---|---|---
RA [rad] | $0.375$ | $3.026$ | $4.154$
DEC [rad] | $-0.363$ | $-0.282$ | $-0.972$
$\Psi$ [rad] | $1.738$ | $4.552$ | $2.071$
$\iota$ [rad] | $1.645$ | $1.143$ | $1.261$
$d_{\rm{L}}$ [Mpc] | $453$ | $453$ | $453$
$z$ | $0.098$ | $0.098$ | $0.098$
$m_{1,\rm src}[M_{\odot}]$ | $10.0$ | $30.0$ | 60.0, 10*
$m_{2,\rm src}[M_{\odot}]$ | $5.0$ | $15.0$ | 30.0, 5*
$t_{\rm c}$ [GPS] | $1120381489.604$ | $1124867056.193$ | $1122192865.797$
phase | $3.472$ | $1.098$ | $3.441$
spin $a_{1}$ | $0.729$ | $0.729$ | $0.729$
spin $a_{2}$ | $0.620$ | $0.620$ | $0.620$

SNR results for these signals are summarized in table 2. The SNRs match
between gwfish and bilby with expected minor deviations, for example, because
gwfish simulates the effects of the rotation of Earth on the waveform while
bilby does not.

Table 2: The comparison between SNR for the four different signals (see Table
1) between bilby and gwfish. The IMRPhenomD waveform model was used with
$f_{\rm{min}}=10$ Hz and $f_{\rm{max}}=1024$ Hz. The SNRs reported here are
the sum (in quadrature) over ET’s three components. SNR
---
Signal | bilby | gwfish
Signal 1 | $72.1$ | $72.1$
Signal 2 | $219$ | $219$
Signal 3 | $1038$ | $1038$
Signal 3* | $236$ | $236$

For the comparison of PEs between bilby and gwfish, we use the same signals of
table 1. Posterior sampling requires to specify a prior on the parameters (at
the very least a uniform prior to constrain the reasonable/physical range of
parameter values). We used:

  * $\bullet$

uniform priors for mass ratio and chirp mass

  * $\bullet$

uniform prior for the polarization angle $\Psi$

  * $\bullet$

uniform prior on the merger rate in comoving volume (fixing the prior on the
luminosity distance $d_{\rm{L}}$)

  * $\bullet$

isotropic distribution in the sky (fixing the priors on right ascension RA and
declination DEC)

  * $\bullet$

uniform prior for the cosine of the inclination angle $\iota$

The other 8 parameters (6 spin parameters, phase and time) are given a delta
prior and therefore are not taken into consideration for posterior sampling
and are excluded also from the Fisher matrices. bilby results are obtained
using the ’PEsummary’ package [40]. In table 3, we report the recovered median
values for the parameters and the relative errors at the $90\%$ confidence
level for bilby. The parameter values reported for gwfish are the true
waveform parameters. gwfish outputs $1\sigma$ PE errors, which are converted
to the $90\%$ confidence level in Table 3, which can be done in Gaussian
approximation by multiplying the errors by $1.645$. Moreover, under the gwfish
results column, for each signal, we have two results per parameter: the upper
one is obtained using $f_{\rm min}=10$ Hz, the lower one with $f_{\rm min}=2$
Hz. Starting from 2 Hz yields smaller uncertainties on the parameters by about
a factor 2.

Generally, there is good agreement between gwfish and bilby PE errors when
taking gwfish results with 10 Hz low-frequency cutoff. For example, mass
errors typically match within a factor 2. There are exceptions, which are all
connected to multi-modal posterior distributions, e.g., several errors
estimated with bilby are significantly larger since they include the
contribution of a second mode, which the Gaussian likelihood approximation
cannot reproduce. Interesting to point out is that at least for these 4
signals, which all have very high SNR, $\iota$, $d_{\rm L}$ are both estimated
well. This tells us that ET, as suspected in the past due to its triangular
configuration, is able to break the inclination-distance degeneracy. It needs
to be noted though that our small and high-SNR sample of signals is not
sufficient to draw general conclusions about ET PE capabilities, and errors
might also increase when spins, merger phase and time are included in the
analysis. What we can conclude from this table is that correlations between
model parameters do not necessarily invalidate the gwfish PE errors, which
means that Fisher-matrix based analyses can produce useful estimates of PE
errors.

Table 3: Comparing measurement uncertainties obtained with bilby and gwfish. We use the signals outlined in Table 1 and IMRPhenomD as waveform model ($f_{\rm{min}}=10$ Hz and $f_{\rm{max}}=1024$ Hz). The errors are at $90\%$ confidence/credibility. The gwfish column presents two results for each parameter: the upper one is obtained starting from $f_{\rm{min}}=10$ Hz, the one below is obtained starting from a minimum frequency of 2 Hz. The parameter values reported for gwfish are the true waveform parameters, while the values reported for bilby are the median values of the posterior distribution. | Signal 1 | Signal 2 | Signal 3 | Signal 3*
---|---|---|---|---
Parameter | gwfish | bilby | gwfish | bilby | gwfish | bilby | gwfish | bilby
RA [rad] |  | $0.37^{+0.3}_{-0.3}$
---
$0.37^{+0.19}_{-0.19}$
$0.5^{+5.1}_{-0.11}$ |  | $3.0^{+0.0046}_{-0.0046}$
---
$3.0^{+0.0031}_{-0.0031}$
$1.2^{+1.9}_{-0.009}$ |  | $4.1^{+1.3}_{-1.3}$
---
$4.1^{+0.84}_{-0.84}$
$4.2^{+0.28}_{-0.11}$ |  | $4.1^{+1.6}_{-1.6}$
---
$4.1^{+0.73}_{-0.73}$
$4.6^{+0.45}_{-0.47}$
DEC [rad] |  | $-0.36^{+0.19}_{-0.19}$
---
$-0.36^{+0.12}_{-0.12}$
$-0.38^{+1.1}_{-0.051}$ |  | $-0.28^{+0.0039}_{-0.0039}$
---
$-0.28^{+0.0031}_{-0.0031}$
$-0.82^{+0.54}_{-0.04}$ |  | $-0.97^{+1.3}_{-1.3}$
---
$-0.97^{+0.81}_{-0.81}$
$-1.0^{+0.32}_{-0.087}$ |  | $-0.97^{+1.6}_{-1.6}$
---
$-0.97^{+0.7}_{-0.7}$
$-0.8^{+0.37}_{-0.29}$
$\Psi$ [rad] |  | $1.7^{+0.3}_{-0.3}$
---
$1.7^{+0.19}_{-0.19}$
$1.8^{+0.46}_{-1.4}$ |  | $4.5^{+0.019}_{-0.019}$
---
$4.5^{+0.012}_{-0.012}$
$1.4^{+2.7}_{-0.44}$ |  | $2.1^{+0.85}_{-0.85}$
---
$2.1^{+0.54}_{-0.54}$
$2.0^{+0.064}_{-0.21}$ |  | $2.1^{+1}_{-1}$
---
$2.1^{+0.47}_{-0.47}$
$4.8^{+0.46}_{-3.3}$
$\iota$ [rad] |  | $1.6^{+0.12}_{-0.12}$
---
$1.6^{+0.076}_{-0.076}$
$1.7^{+0.047}_{-0.046}$ |  | $1.1^{+0.045}_{-0.045}$
---
$1.1^{+0.031}_{-0.031}$
$1.1^{+0.048}_{-0.06}$ |  | $1.3^{+0.0089}_{-0.0089}$
---
$1.3^{+0.0057}_{-0.0057}$
$1.3^{+0.001}_{-0.002}$ |  | $1.3^{+0.012}_{-0.012}$
---
$1.3^{+0.0071}_{-0.0071}$
$1.3^{+0.006}_{-0.007}$
$d_{\rm{L}}$ [Mpc] |  | $453^{+400}_{-400}$
---
$453^{+246}_{-246}$
$576^{+66}_{-120}$ |  | $453^{+43}_{-43}$
---
$453^{+29}_{-29}$
$477^{+53}_{-44}$ |  | $453^{+8}_{-8}$
---
$453^{+5}_{-5}$
$453^{+1.1}_{-1.2}$ |  | $453^{+10}_{-10}$
---
$453^{+5}_{-5}$
$453^{+4.7}_{-4.6}$
$\text{m}_{1,\rm src}[M_{\odot}]$ |  | $10^{+0.095}_{-0.095}$
---
$10^{+0.036}_{-0.036}$
$9.4^{+0.74}_{-2.5}$ |  | $30^{+0.22}_{-0.22}$
---
$30^{+0.074}_{-0.074}$
$30^{+0.62}_{-0.63}$ |  | $60^{+0.13}_{-0.13}$
---
$60^{+0.058}_{-0.058}$
$60^{+0.095}_{-0.092}$ |  | $10^{+0.11}_{-0.11}$
---
$10^{+0.061}_{-0.061}$
$10^{+0.11}_{-0.1}$
$\text{m}_{2,\rm src}[M_{\odot}]$ |  | $5.0^{+0.042}_{-0.042}$
---
$5.0^{+0.016}_{-0.016}$
$5.2^{+1.7}_{-0.36}$ |  | $15^{+0.099}_{-0.099}$
---
$15^{+0.033}_{-0.033}$
$15^{+0.17}_{-0.17}$ |  | $30^{+0.06}_{-0.06}$
---
$30^{+0.026}_{-0.026}$
$30^{+0.04}_{-0.04}$ |  | $5.0^{+0.046}_{-0.046}$
---
$5.0^{+0.026}_{-0.026}$
$5.0^{+0.046}_{-0.045}$

Corner plots of the posterior samples with bilby are shown in the appendix in
Figures 6 to 9.

###  III.2 TaylorF2 versus IMRPhenomD versus IMRPhenomXPHM

Next, we compare results obtained with two approximants: TaylorF2 (see [41]
and [42]) and IMRPhenomD (see [43] and [44]). The TaylorF2 model is the
simplest frequency-domain approximant, which describes the inspiral phase, and
accurately only the early inspiral, which is sufficient for estimates, e.g.,
of the chirp mass and BNS localization through Earth’s rotation. We apply a
high-frequency cutoff at 4 times the frequency of the innermost stable
circular orbit $f_{\rm isco}=6^{-1.5}\pi^{-1}c^{3}/(GM)$ to approximately
match the SNRs obtained with other waveform models. The IMRPhenomD model
extends over the full inspiral-merger-ringdown of a compact-binary signal, but
there are still constraints on the physics described by this model. For
example, it only models the fundamental mode of the GW signal and systems with
aligned spins. The purpose of this comparison is to show how a waveform model
affects parameter estimation, even at the simple level of Fisher-matrix based
error estimates.

Figure 2: Comparing TaylorF2 and IMRPhenomD waveforms using ET and Signal 2.
Normalized strain in frequency domain for both the waveform models is plotted
and compared to the used detector sensitivity. A cut at $4f_{\rm{isco}}$ is
applied to the TaylorF2 waveform.

To highlight the model differences we display in Figure 2 the normalized
amplitudes for the two waveforms for Signal 2 (see Table 1) adjusted to units
Hz-1/2 in comparison also with the ET sensitivity curve. It can be observed
that the IMRPhenomD model shows the characteristic bump towards the merger
frequency and provides a more natural cut off at high frequency.

To compare the two approximants, we start with 1000 BBH signals and 1000 BNS
signals. We constrain the BNS redshifts to a maximum of 0.4 so that the SNR
distributions for BBHs and BNSs are similar. From these, 927 BBH and 966 BNS
signals are seen with a minimum SNR of 8.0 using TaylorF2, IMRPhenomD, or
IMRPhenomXPHM LALSimulation approximants. The distribution of SNR values of
detected signals goes up to 200 peaking at 20. The signals are analyzed with a
frequency resolution of 1/32 Hz. The analysis is carried out with ET alone
over a frequency range from 2 Hz to 1024 Hz. The results are shown in Figure
3. We plot the ratios of the PE errors recovered for each of the parameters
with TaylorF2 and IMRPhenomD as histograms over the detected signals. The BBH
results are in gray, while the BNS ones are plotted in red. The upper graph
shows the ratio of SNRs, while the other histograms show the ratio of the
uncertainties on the other parameters. For BBH systems, the PE errors of the
individual masses are underestimated with TaylorF2 (assuming that IMRPhenomD
generally yields more accurate results), while all other parameters considered
here are weakly biased towards larger errors compared to IMRPhenomD. For BNS
signals there is no significant difference between the two considered
approximants: the SNRs are coincident, and only the sky-localisation errors
are slightly underestimated with TaylorF2.

Figure 3: Comparison between TaylorF2 and IMRPhenomD approximants for BBH
signals (in gray) and for BNS signals (in red). We analyze 927 BBH signals and
966 BNS systems observed by ET with SNR $>$ 8\. The frequency range is from 2
Hz to 1024 Hz. The frequency resolution used is 1/32 Hz. We applied a cut at
$4f_{\rm{isco}}$ for the TaylorF2 waveform. The upper histogram shows the
ratio between the SNRs. The histograms below show the ratio of the PE errors.
The vertical dashed black line marks the equal ratio case. Figure 4:
Comparison between IMRPhenomXPHM and IMRPhenomD approximants for BBH signals
(in gray) and for BNS signals (in red). Settings and number of signals are the
same as in Figure 3.

We show a similar analysis for the IMRPhenomXPHM and IMRPhenomD approximants.
The IMRPhenomXPHM approximant can describe precession (not considered in this
study) and higher-order modes. The latter is well-known to have important
consequences for BBH analyses (see, e.g., [45, 46]). Here, a significant
improvement of PE can also be seen for a significant fraction of BNS signals.
For about 5% of BNS, the sky-localization error is significantly reduced by up
to a factor 5 using IMRPhenomXPHM. For about 10% of BNS, the distance error is
significantly reduced by up to a few orders of magnitude. What this means is
that some BNS signals show strong correlations in the likelihood between
distance and inclination angle using IMRPhenomD, and the higher-order modes
described with IMRPhenomXPHM lead to a suppression of these correlations.

##  IV Einstein Telescope

The Einstein Telescope is a proposed underground infrastructure in Europe to
host future-generations of GW detectors. With respect to its PE capabilities,
three features stand out, which will make it a unique facility among future
detectors:

  * $\bullet$

It has the geometry of an equilateral triangle and laser-interferometer arms
launching from each of its vertices. It is expected that this triangular
configuration will give it unique capabilities concerning the measurement of
GW polarizations and parameters correlated with GW polarizations.

  * $\bullet$

It will have 10 km arm length, which, together with several technology
improvements with respect to current detectors, leads to a greatly improved
peak sensitivity by about a factor 10. Binaries of compact objects will be
detected along the cosmic history; BNS and BBH mergers are both expected to be
detected with a rate of $10^{5}$ per year. The enormous increase of SNR for
close events will enable accurate parameter estimation, enabling great
precision in GW astronomy.

  * $\bullet$

It will consist of a total of 3 pairs of laser interferometers. Each pair
includes a low-frequency interferometer and a high-frequency interferometer
together covering an observation band from 3 Hz to several kHz. The
observation is thus extended to much lower frequencies compared to current GW
detectors, which leads to longer observation times of individual compact-
binary signals, and it will also make it possible to observe larger-mass
compact binaries well into the regime of intermediate-mass black-hole
binaries. Long observation times enable stand-alone sky-location measurements
by ET of neutron-star binaries exploiting amplitude modulations due to Earth’s
rotation. Due to the fact that the observed mass is the redshifted source-
frame mass, i.e., multiplication by $(1+z)$, the low-frequency sensitivity
will make it possible to access BBH mergers in the early Universe where
primordial BHs can be discriminated from BHs with stellar origin [e.g. 47,
48].

In the following, we present a few reference PE results obtained from gwfish,
which can be used for comparison with other simulation software. We focus on
sky-localization results, since these are more difficult to simulate
accurately. A single detector can accurately localize a GW source if the
observation is long enough so that the motion of the detector leads to
amplitude and phase modulations of GW signals, which depend on the direction
of the source. This is the case for BNSs, which are observed for several hours
with ET. Even lighter NSBHs and BBHs can be localized this way by ET. For
specific signals as well as for the average localization by Earth rotation,
the localization errors depend on the latitude of the GW detector. Sky-
localization is greatly improved if additional detectors are present to form a
network. In the latter case, the information comes from measurements of phase
differences between signals observed in different detectors (standard
triangulation technique) [49].

We show in Figure 5 the sky-localization results for BNSs as $90\%$ credible
region. The results are presented as fraction of signals localized within a
certain area in the sky. BNSs are considered at various fixed luminosity
distances; 1000 BNS per curve. The BNS masses are set to 1.4+1.4$M_{\odot}$.
The plots confirm that ET has unique sky-localization capabilities for BNSs,
which would enable multi-messenger studies without the need of other
detectors, but it is also clear that the localization capabilities of a
network are far superior even if only one other detector like Cosmic Explorer
(CE) is added. Nearer systems are better localized and the network of ET+2CE
manages to localize almost all of the injected signals within 10 deg2. We find
our results to be consistent with previously published results, e.g., [15,
50].

Figure 5: Sky-localization capabilities (90% credible regions) of future
detector networks including ET and CE. We used the LIGO sites as locations for
the two CE detectors. Left: ET alone, middle: ET and one CE detector at the
LIGO Hanford site, right: ET with two CE detectors at the two LIGO sites. The
sensitivity models used here is the ET sensitivity formerly known as ET-D
[51], and the 40 km CE sensitivity [9].

Next, we present an analysis of the 85 signals of GWTC-3 [52, 53, 54]. Only
the merger time, masses and distance were taken from the catalogue. Angular
parameters were drawn randomly, spins were set to zero. Observation of these
signals with ET was simulated. The simulation was done for 10 realizations of
the angular waveform parameters for each signal. Table 4 provides a summary of
the median SNRs and sky-localization that would have been achieved with ET.
The table only contains the signals whose sky-localization was within 1000
deg2 for all random realizations of the angular variables of the waveform. As
expected, these signals are among the lowest-mass systems in GWTC-3, which
means that ET would have observed them for a long time.

| GW170817 | GW190425 | GW190707
---|---|---|---
Median SNR | 710 | 190 | 230
Median sky-location error [deg2] (90% conf.) | 4.0 | 67 | 34
Masses [$M_{\odot}$] | 1.46, 1.27 | 2, 1.4 | 12, 8.4
Luminosity distance [Mpc] | 40 | 160 | 800
Time until merger starting from 3 Hz [s] | $2.6\cdot 10^{4}$ | $1.8\cdot 10^{4}$ | 930
Table 4: Median SNRs and sky-localization errors for simulated observations of
GWTC-3 signals with ET. Only masses, merger times, and distances were taken
from GWTC-3. The sky-localization errors were evaluated for several random
realizations, and only the signals with error $<1000$ deg2 for all
realizations are shown in the table. The waveform used here is IMRPhenomXPHM
to include the potentially important higher-order modes, but for the signals
reported in the table, the sky-localization results would have been the same
using TaylorF2 or IMRPhenomD. The localization even of larger BBHs can be
$<100$ deg2, but it depends strongly on the values of angular parameters.

##  V Conclusions

In this paper, we have presented gwfish; a new Fisher-matrix code to
investigate parameter-estimation capabilities of future GW detectors. gwfish
stands out as a robust and flexible detector simulation software:

  * $\bullet$

It implements modeling of the so-called $A,E,T$ channels of LISA-type space-
borne detectors, of LGWA, and of triangular and L-shape terrestrial GW
detectors;

  * $\bullet$

It can use all LALSimulation approximants supporting precession as well as
tidal deformation;

  * $\bullet$

It incorporates the motion of detectors, e.g., because of Earth/Moon rotation,
or the orbits of satellites around the Sun, and it simulates the corresponding
amplitude and phase modulations of GW signals;

  * $\bullet$

It can simulate multi-band observations;

  * $\bullet$

It takes several measures to reduce the impact of numerical errors on Fisher-
matrix PE errors.

Certain analyses are currently not possible with state-of-the-art detector
simulation and Bayesian analysis software like bilby due to practical
constraints on available computational resources. For example, if hundreds of
thousands of GW signals need to be analyzed, sometimes repeatedly to study the
effect of changes in population or waveform models, then Fisher-matrix codes
like gwfish are the best option. Important is to understand the limitations of
Gaussian approximations of the likelihood. We have shown a few comparisons in
this article between bilby and gwfish, and we generally find a good match
between the two with exceptions that can be explained. In the cases studied
here, discrepancies emerged when multiple modes of the posterior distribution
contributed significantly to the PE errors. Tighter priors would significantly
reduce the mismatch.

We also investigated the impact of the choice of waveform model, i.e.,
TaylorF2 vs IMRPhenomD vs IMRPhenomXPHM. For BBH signals, all three
approximants produce significantly different results. For the BNS analyses
presented here, namely SNR and sky-localization, TaylorF2 vs IMRPhenomD
produce almost identical results. However, for a significant fraction of BNS
signals, PE errors with IMRPhenomXPHM are significantly reduced compared to
IMRPhenomD. The reduction is strongest for distance errors (10% of signals
significantly affected), because distance errors can be strongly correlated
with errors of the inclination angle. The IMRPhenomXPHM approximant takes
significantly longer to evaluate than TaylorF2 and IMRPhenomD (the latter two
with similar evaluation times) when using the LALSimulation package with
Python/SWIG [55]. We conclude that the choice of approximant is generally more
important for BBH analyses. For BNS signals, broad population studies can
still be carried out with IMRPhenomD or TaylorF2, but results especially for
distance and inclination errors are more reliable for a significant fraction
of signals when using IMRPhenomXPHM. Needless to say, it always depends on the
type of analysis one wants to perform whether a certain waveform model is
adequate or not. Considering precession and higher-order modes might be
negligible effects in some analyses, and crucial in others. The gwfish user
must be familiar with the characteristics of the standard waveform models.

The scientific results presented in this paper confirm the long-standing claim
that the triangular configuration of ET makes it possible to measure GW
polarizations and thereby break certain model degeneracies. For example, for
high-SNR BBH signals, we found that the inclination angle is measured
precisely with ET, which also improves the estimates of source distance. We
leave it for a future analysis whether this ability extends to lower-SNR
signals, which could be of great benefit for cosmological studies. We also
found that ET itself would have estimated the sky location of the low-mass
GWTC-3 signals to within 100 deg2 at 90% confidence (depending on the exact
sky location, the localization of GW170817 can be within 1 deg2). This shows
that multi-messenger astronomy with ET as sole GW detector is possible, but
unlocking the full potential of multi-messenger astronomy with GWs would
require at least another GW detector with a similar sensitivity.

## Acknowledgments

We thank Michele Mancarella, Francesco Iacovelli, Stefano Foffa, and Michele
Maggiore for discussions about the Fisher-matrix technique and for comparisons
of results with their simulation code. We thank Stefano Bagnasco and the IT
team of University of Turin for support with computational resources. BG is
supported by the Italian Ministry of Education, University and Research within
the PRIN 2017 Research Program Framework, n. 2017SYRTCN. MB, AM and SR
acknowledge financial support from MIUR (PRIN 2020 grant 2020KB33TP_001). BB
acknowledges financial support from MIUR (PRIN 2017 grant 20179ZF5KS).

## References

  * Maggiore et al. [2020] M. Maggiore, C. V. D. Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse, S. Foffa, J. García-Bellido, et al., Journal of Cosmology and Astroparticle Physics 2020, 050 (2020), URL https://doi.org/10.1088%2F1475-7516%2F2020%2F03%2F050.
  * Kawamura et al. [2021] S. Kawamura, M. Ando, N. Seto, S. Sato, M. Musha, I. Kawano, J. Yokoyama, T. Tanaka, K. Ioka, T. Akutsu, et al., Progress of Theoretical and Experimental Physics 2021 (2021), ISSN 2050-3911, URL https://doi.org/10.1093/ptep/ptab019.
  * Acernese et al. [2014] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, et al., Classical and Quantum Gravity 32, 024001 (2014), URL https://doi.org/10.1088%2F0264-9381%2F32%2F2%2F024001.
  * Aasi et al. [2015] J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al., Classical and Quantum Gravity 32, 074001 (2015), URL https://doi.org/10.1088%2F0264-9381%2F32%2F7%2F074001.
  * Souradeep [2016] T. Souradeep, Resonance 21, 225 (2016), URL https://doi.org/10.1007/s12045-016-0316-6.
  * Akutsu et al. [2019] T. Akutsu, M. Ando, K. Arai, Y. Arai, S. Araki, A. Araya, N. Aritomi, H. Asada, Y. Aso, S. Atsuta, et al., Nature Astronomy 3, 35 (2019), ISSN 2397-3366, URL https://www.nature.com/articles/s41550-018-0658-y.
  * Punturo et al. [2010] M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia, M. Beker, et al., Classical and Quantum Gravity 27, 194002 (2010), URL http://stacks.iop.org/0264-9381/27/i=19/a=194002.
  * ET Steering Committee [2020] ET Steering Committee, available from European Gravitational Observatory, document number ET-0007B-20 (2020), URL https://apps.et-gw.eu/tds/ql/?c=15418.
  * Evans et al. [2021] M. Evans, R. X. Adhikari, C. Afle, S. W. Ballmer, S. Biscoveanu, S. Borhanian, D. A. Brown, Y. Chen, R. Eisenstein, A. Gruson, et al., _A horizon study for cosmic explorer: Science, observatories, and community_ (2021), URL https://arxiv.org/abs/2109.09882.
  * Amaro-Seoane et al. [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al., arXiv e-prints (2017), URL https://arxiv.org/abs/1702.00786.
  * Hall et al. [2021] E. D. Hall, K. Kuns, J. R. Smith, Y. Bai, C. Wipf, S. Biscans, R. X. Adhikari, K. Arai, S. Ballmer, L. Barsotti, et al., Phys. Rev. D 103, 122004 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.122004.
  * Sesana [2016] A. Sesana, Phys. Rev. Lett. 116, 231102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.231102.
  * Vitale [2016] S. Vitale, Phys. Rev. Lett. 117, 051102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.051102.
  * Grimm and Harms [2020] S. Grimm and J. Harms, Phys. Rev. D 102, 022007 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.022007.
  * Chan et al. [2018] M. L. Chan, C. Messenger, I. S. Heng, and M. Hendry, Phys. Rev. D 97, 123014 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.97.123014.
  * Borhanian [2021] S. Borhanian, Classical and Quantum Gravity 38, 175014 (2021), URL https://doi.org/10.1088/1361-6382/ac1618.
  * Cutler [1998] C. Cutler, Phys. Rev. D 57, 7089 (1998), URL https://link.aps.org/doi/10.1103/PhysRevD.57.7089.
  * Nitz and Canton [2021] A. H. Nitz and T. D. Canton, The Astrophysical Journal Letters 917, L27 (2021), URL https://doi.org/10.3847/2041-8213/ac1a75.
  * Vallisneri [2008] M. Vallisneri, Phys. Rev. D 77, 042001 (2008), URL https://link.aps.org/doi/10.1103/PhysRevD.77.042001.
  * Rodriguez et al. [2013] C. L. Rodriguez, B. Farr, W. M. Farr, and I. Mandel, Phys. Rev. D 88, 084013 (2013), URL https://link.aps.org/doi/10.1103/PhysRevD.88.084013.
  * Armstrong et al. [1999] J. W. Armstrong, F. B. Estabrook, and M. Tinto, The Astrophysical Journal 527, 814 (1999), URL https://doi.org/10.1086/308110.
  * Prince et al. [2002] T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, Phys. Rev. D 66, 122002 (2002), URL https://link.aps.org/doi/10.1103/PhysRevD.66.122002.
  * Vallisneri [2005] M. Vallisneri, Phys. Rev. D 71, 022001 (2005), URL https://link.aps.org/doi/10.1103/PhysRevD.71.022001.
  * Hinderer et al. [2010] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Phys. Rev. D 81, 123016 (2010), URL http://link.aps.org/doi/10.1103/PhysRevD.81.123016.
  * Ronchini et al. [2022] S. Ronchini, M. Branchesi, G. Oganesyan, B. Banerjee, U. Dupletsa, G. Ghirlanda, J. Harms, M. Mapelli, and F. Santoliquido, arXiv e-prints arXiv:2204.01746 (2022), eprint 2204.01746.
  * Goncharov et al. [2022] B. Goncharov, A. H. Nitz, and J. Harms, _Utilizing the null stream of einstein telescope_ (2022), URL https://arxiv.org/abs/2204.08533.
  * Schilling [1997] R. Schilling, Classical and Quantum Gravity 14, 1513 (1997), URL https://doi.org/10.1088/0264-9381/14/6/020.
  * Abbott et al. [2017] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, K. Ackley, C. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, C. Affeldt, et al., Classical and Quantum Gravity 34, 044001 (2017), URL http://stacks.iop.org/0264-9381/34/i=4/a=044001.
  * Essick et al. [2017] R. Essick, S. Vitale, and M. Evans, Phys. Rev. D 96, 084004 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.084004.
  * Droz et al. [1999] S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, Phys. Rev. D 59, 124016 (1999), URL https://link.aps.org/doi/10.1103/PhysRevD.59.124016.
  * Harms et al. [2013] J. Harms, B. J. J. Slagmolen, R. X. Adhikari, M. C. Miller, M. Evans, Y. Chen, H. Müller, and M. Ando, Phys. Rev. D 88, 122003 (2013), URL http://link.aps.org/doi/10.1103/PhysRevD.88.122003.
  * Sathyaprakash and Schutz [2009] B. S. Sathyaprakash and B. F. Schutz, Living Reviews in Relativity 12, 2 (2009), ISSN 1433-8351, URL https://doi.org/10.12942/lrr-2009-2.
  * Gürsel and Tinto [1989] Y. Gürsel and M. Tinto, Phys. Rev. D 40, 3884 (1989), URL https://link.aps.org/doi/10.1103/PhysRevD.40.3884.
  * Regimbau et al. [2012] T. Regimbau, T. Dent, W. Del Pozzo, S. Giampanis, T. G. F. Li, C. Robinson, C. Van Den Broeck, D. Meacher, C. Rodriguez, B. S. Sathyaprakash, et al., Phys. Rev. D 86, 122001 (2012), URL https://link.aps.org/doi/10.1103/PhysRevD.86.122001.
  * Parida et al. [2016] A. Parida, S. Mitra, and S. Jhingan, Journal of Cosmology and Astroparticle Physics 2016, 024 (2016), URL https://doi.org/10.1088/1475-7516/2016/04/024.
  * Press et al. [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge University Press, 2007).
  * LIGO Scientific Collaboration [2018] LIGO Scientific Collaboration, _LIGO Algorithm Library - LALSuite_ , free software (GPL) (2018).
  * Ashton et al. [2019] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, et al., The Astrophysical Journal Supplement Series 241, 27 (2019), URL https://doi.org/10.3847%2F1538-4365%2Fab06fc.
  * Vallisneri [2011] M. Vallisneri, Phys. Rev. Lett. 107, 191104 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.191104.
  * Hoy and Raymond [2021] C. Hoy and V. Raymond, SoftwareX 15, 100765 (2021), eprint 2006.06639.
  * Buonanno et al. [2009] A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Physical Review D 80 (2009), ISSN 1550-2368, URL http://dx.doi.org/10.1103/PhysRevD.80.084043.
  * Isoyama et al. [2021] S. Isoyama, R. Sturani, and H. Nakano, Handbook of Gravitational Wave Astronomy p. 1–49 (2021), URL http://dx.doi.org/10.1007/978-981-15-4702-7_31-1.
  * Husa et al. [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. J. Forteza, and A. Bohé, Physical Review D 93 (2016), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.93.044006.
  * Khan et al. [2016] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. J. Forteza, and A. Bohé, Physical Review D 93 (2016), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.93.044007.
  * Broeck and Sengupta [2007] C. V. D. Broeck and A. S. Sengupta, Classical and Quantum Gravity 24, 1089 (2007), URL https://doi.org/10.1088%2F0264-9381%2F24%2F5%2F005.
  * Khan et al. [2020] S. Khan, F. Ohme, K. Chatziioannou, and M. Hannam, Phys. Rev. D 101, 024056 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.024056.
  * Ng et al. [2021] K. K. Y. Ng, S. Chen, B. Goncharov, U. Dupletsa, S. Borhanian, M. Branchesi, J. Harms, M. Maggiore, B. S. Sathyaprakash, and S. Vitale, arXiv e-prints arXiv:2108.07276 (2021), eprint 2108.07276.
  * Luca et al. [2021] V. D. Luca, G. Franciolini, P. Pani, and A. Riotto, Journal of Cosmology and Astroparticle Physics 2021, 003 (2021), URL https://doi.org/10.1088/1475-7516/2021/05/003.
  * Wen and Chen [2010] L. Wen and Y. Chen, Phys. Rev. D 81, 082001 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.81.082001.
  * Li et al. [2022] Y. Li, I. S. Heng, M. L. Chan, C. Messenger, and X. Fan, Physical Review D 105 (2022), URL https://doi.org/10.1103%2Fphysrevd.105.043010.
  * Hild et al. [2011] S. Hild, M. Abernathy, F. Acernese, P. Amaro-Seoane, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia, M. Beker, et al., Classical and Quantum Gravity 28, 094013 (2011), URL http://stacks.iop.org/0264-9381/28/i=9/a=094013.
  * Abbott et al. [2019] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 9, 031040 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.031040.
  * Abbott et al. [2021] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 11, 021053 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.021053.
  * The LIGO Scientific Collaboration et al. [2021] The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. X. Adhikari, et al., _GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run_ (2021), URL https://arxiv.org/abs/2111.03606.
  * Wette [2020] K. Wette, SoftwareX 12, 100634 (2020), ISSN 2352-7110, URL https://www.sciencedirect.com/science/article/pii/S2352711020303472.

##  Appendix A Posterior samples

Figures 6, 7, 8 and 9 show corner plots of the posteriors obtained with bilby
for signals 1, 2, 3 and 3* respectively, whose parameters are described in
table 1.

Figure 6: Posterior sampling output for Signal 1 ($5M_{\odot}+10M_{\odot}$).
Figure 7: Posterior sampling output for Signal 2 ($15M_{\odot}+30M_{\odot}$).
Figure 8: Posterior sampling output for Signal 3 ($30M_{\odot}+60M_{\odot}$).
Figure 9: Posterior sampling output for Signal 3* ($5M_{\odot}+10M_{\odot}$).


# Semicocycle discontinuities for substitutions and reverse-reading automata

Gandhar Joshi  School of Mathematics and Statistics, The Open University, U.K.
<EMAIL_ADDRESS>and  Reem Yassawi  School of Mathematical Sciences,
Queen Mary University of London, UK<EMAIL_ADDRESS>
###### Abstract.

In this article we define the semigroup associated to a substitution. We use
it to construct a minimal automaton which generates a substitution sequence u
in reverse reading. We show, in the case where the substitution has a
coincidence, that this automaton completely describes the semicocycle
discontinuities of u.

This work was supported by the EPSRC grant number EP/V007459/1.

In Memoriam: Uwe Grimm died unexpectedly in October 2021. He was the first
author’s supervisor and the second author’s close colleague and dear friend.
His loss is deeply felt by both; we dedicate this paper to his memory.

##  1\. Introduction

Automatic sequences are codings of constant length-$\ell$ substitutional fixed
points. The latter generate discrete dynamical systems whose study is
ubiquitous in the literature. Cobham [2] showed that an equivalent definition
of an automatic sequence is that it is a sequence generated by a finite
automaton with output. To generate $u_{n}$, the $n$-th entry of the sequence,
one feeds the base-$\ell$ expansion $(n)_{\ell}$ of $n$ into the automaton, to
arrive at a final state, and $u_{n}$ is the coding (“output”) of the final
state. One has to set the order of reading of $(n)_{\ell}$, i.e., whether one
reads starting with the most significant digit, in which case the reading is
called direct, or starting with the least significant digit, in which case the
reading is reverse. A fixed automatic sequence can be generated in either
direct or reverse reading, but the generating automata will be different.

The relevance of the (minimal) direct-reading automaton in dynamics has long
been understood. Indeed, Cobham’s proof tells us that we can define this
automaton directly in terms of the substitution rule. Furthermore, this
direct-reading automaton also enables us to locate the irregular fibres of the
maximal equicontinuous factor [3], and this is useful as these fibres often
drive the interesting dynamics in the shift. Cobham defined automatic
sequences in terms of direct reading automata, and he did not consider
reverse-reading automata. Eilenberg identified that reverse-reading automata
gave information about the $\ell$-kernel of a sequence u, that is, the set of
all subsequences of u whose indices run along an arithmetic progression and
such that its difference is a power of $\ell$. Furthermore, the cardinality of
the $\ell$-kernel is also the size of a minimal automaton that generates u in
reverse reading. The reverse-reading automaton has been used extensively to
generate algebraic characterisations of automatic sequences, in particular
when their entries belong to a field or even a ring. However there has been
little or no use of the reverse-reading automaton in dynamics.

This paper arose out of a desire to understand the dynamical meaning of the
$\ell$-kernel, and also to understand what the reverse-reading automaton tells
us about the dynamical system generated by u. It turns out that for a special
class of shift-dynamical systems, namely the Toeplitz shifts, the
$\ell$-kernel is relevant. These shifts are generated as the shift-orbit
closure of a Toeplitz sequence, i.e., one which is constructed by filling in
entries of one arithmetic progression at a time, with a constant symbol, in a
way such that eventually the sequence is completely specified. The first such
sequences were studied by Oxtoby [9], who defined a Toeplitz sequence that
generates a minimal shift supporting two invariant measures.

At the other end, Jacobs and Keane [6], introduced what is possibly the
simplest nontrivial Toeplitz sequence, the well-known period-doubling
sequence. It is a substitutional fixed point, and generates a uniquely ergodic
shift which has very little “irregularity”, as follows. Let $(X,\sigma)$ be a
Toeplitz shift. Not all elements in $X$ are Toeplitz sequences. How far away a
sequence $(x_{n})_{n\in\mathbb{I}}$ is from being Toeplitz is a function of
the “size” of the discontinuities of the semicocycle that it defines. Here the
semicocycle is the map $n\in\mathbb{I}\mapsto x_{n}\in\mathcal{A}$, where the
topology on $\mathbb{I}$ is inherited by the topology of the maximal
equicontinuous factor of $(X,\sigma)$, and in this paper, $\mathcal{A}$ is a
finite alphabet with the discrete topology. There is a discontinuity each time
an index does not belong to an arithmetic progression that has been filled in
with a constant letter in $\mathcal{A}$. The period-doubling sequence defines
a semicocycle on $\mathbb{Z}$ with exactly one discontinuity.

Williams [10] was the first person to study Toeplitz shifts systematically.
Given a point $\textit{{x}}=(x_{n})$ in $X$, she partitions its indices up
into two parts, the periodic part $\mbox{Per}_{k}(\textit{{x}})$, and the
aperiodic part $\mbox{Aper}(\textit{{x}})$. The periodic part is the set of
indices that have been filled using arithmetic progressions, the aperiodic
part is what is left over. Thus, if x is actually a Toeplitz sequence, then
$\mbox{Aper}(\textit{{x}})$ is empty, and the semicocycle is continuous. The
existence of a large $\mbox{Aper}(\textit{{x}})$ is what makes the dynamics
interesting, e.g., it can lead to the existence of several invariant measures,
measurable but not continuous eigenvalues, positive entropy, and so on.

Our main result is Theorem 18. It tells us given an automatic sequence u which
belongs to a Toeplitz shift, the reverse-reading automaton which generates u
will give us complete information about the indices which belong to
$\mbox{Aper}(\textit{{u}})$, and hence the set of discontinuities of the
semicocycle u. To prove Theorem 18, we first construct a suitable automaton
that generates u in reverse reading, in Theorem 12. This automaton is minimal,
so it can also be obtained using Eilenberg’s classic construction, but we
label the states in a different manner using the semigroup of the
corresponding substitution (Definition 9). Finally, consider a formal power
series $f(x)$ whose coefficients form an automatic sequence. Reverse-reading
automata are instrumental in expressing $f(x)$ as either a root of a
polynomial, as given by Christol’s theorem, or as a diagonal of a higher
dimensional rational function, as given by Furstenberg’s theorem; for an
exposition of either result, see [1]. In future work, we will use the reverse-
reading automaton that we have constructed to translate certain dynamical
features of substitution shifts to algebraic properties of their annihilating
polynomials.

##  2\. Preliminaries

Let $\mathbb{N}$ and $\mathbb{N}_{0}$ denote the positive integers and the
non-negative integers respectively. Given $n\in\mathbb{N}_{0}$, let
$(n)_{\ell}$ denote the base-$\ell$ expansion of $n$, where we drop all
leading zeros. If $n\in-\mathbb{N}$, its base-$\ell$ expansion is an infinite
sequence which is ultimately constant with entries equal to $\ell-1$.

Let $\mathcal{A}$ be a finite alphabet. In this article, the indexing set
$\mathbb{I}$ equals $\mathbb{N}_{0}$ or $\mathbb{Z}$. An element in
$\mathcal{A}^{\mathbb{I}}$ is written
$\textit{{u}}=\left(u_{n}\right)_{n\in\mathbb{I}}$. The (left) shift map
$\sigma:\mathcal{A}^{\mathbb{I}}\rightarrow\mathcal{A}^{\mathbb{I}}$ is the
map defined as $(\sigma(\textit{{u}}))_{n}:=u_{n+1}$. Let $\mathcal{A}$ be
endowed with the discrete topology and $\mathcal{A}^{\mathbb{I}}$ with the
product topology; then $\mathcal{A}^{\mathbb{I}}$ is a metrisable Cantor
space. We use the notation $\overline{a}$ to denote a constant sequence whose
entries consist entirely of the letter $a$.

Next we define some finite-state automata with which we will work. We gently
modify existing definitions to allow the generation of two-sided sequences
(see Definition 2).

######  Definition 1.

[Finite automaton] A deterministic finite automaton  (DFA), is a $4$-tuple
$\mathcal{M}=(\mathcal{S},\Sigma,\delta,\mathcal{S}_{0})$, where $\mathcal{S}$
is a finite set (of states), $\mathcal{S}_{0}$ is a set of initial states,
$\Sigma$ is a finite alphabet (the input alphabet), and
$\delta:\mathcal{S}\times\Sigma\rightarrow\mathcal{S}$ (the transition
function).

A deterministic finite automaton with output (DFAO), is a $6$-tuple
$\mathcal{M}=(\mathcal{S},\Sigma,\delta,\mathcal{S}_{0},\mathcal{A},\Omega_{0})$,
where $(\mathcal{S},\Sigma,\delta,\mathcal{S}_{0})$ is a DFA, $\mathcal{A}$ is
a finite alphabet (the output alphabet), and where for each initial state
$s_{0}\in\mathcal{S}_{0}$ there is associated a map $\omega_{0}\in\Omega_{0}$,
$\omega_{0}:\mathcal{S}\rightarrow\mathcal{A}$ (the output function).

The function $\delta$ extends in a natural way to the domain
$\mathcal{S}\times\Sigma^{*}$, where $\Sigma^{*}$ is the set of all finite
words on the alphabet $\Sigma$. One way to do this is to define
$\delta(s,n_{k}\cdots n_{1}n_{0}):=\delta(\delta(s,n_{0}),n_{k}\cdots n_{1})$
recursively. Here, the way we have written things, we first feed $n_{0}$, then
$n_{1}$, etc. To emphasise this, we will write $\delta(s,n_{k}\cdots
n_{1}n_{0})=\delta(s,\stackrel{{\scriptstyle\xleftarrow{\;\;\;\,\quad\qquad}}}{{n_{k}\cdots
n_{1}n_{0}}})$. However, we can also extend $\delta$ to words by starting with
$n_{k}$, and ending with $n_{0}$. In other words we can also define
$\delta(s,n_{k}\cdots n_{1}n_{0}):=\delta(\delta(s,n_{k}),n_{k-1}\cdots
n_{0})$. In this case we will write $\delta(s,n_{k}\cdots
n_{1}n_{0})=\delta(s,\stackrel{{\scriptstyle\xrightarrow{\;\;\;\,\quad\qquad}}}{{n_{k}\cdots
n_{1}n_{0}}})$.

If $n\in\mathbb{N}$, recall that the canonical base-$\ell$ expansion of $n$ is
$(n)_{\ell}=n_{j}\cdots n_{0}$, with $n_{j}\neq 0$. If $n\in-\mathbb{N}$, its
expansion $(n)_{\ell}=\overline{\ell-1}\,n_{j}\cdots n_{0}$ must eventually
equal $\overline{\ell-1}$. For $n\in-\mathbb{N}$, we write its expansion as
$(n)_{\ell}=\overline{\ell-1}\,n_{j}\cdots n_{0}$, where $n_{j}\neq\ell-1$, or
$n_{j}\cdots n_{0}=\emptyset$. We call the finite word $\ell-1\,n_{j}\cdots
n_{0}$ the canonical base-$\ell$ expansion of $n\in-\mathbb{N}$. This is
because we will feed the transition function $\delta$ with base-$\ell$
expansions of integers, and the function $\delta$ only accepts finite words
from $\Sigma^{*}$. Furthermore, we will arrange it so that the canonical
expansion of a negative integer is sufficient for our needs.

In what follows, we will feed the automaton with $(n_{k}\cdots n_{1}n_{0})$
the canonical base-$\ell$ expansion of $n\in\mathbb{I}$, where the digit
$n_{0}$ is the least significant digit and $n_{k}$ is the most significant
digit.

######  Definition 2 (Automaticity).

A sequence $\left(u_{n}\right)_{n\geq 0}$ of elements in $\mathcal{A}$ is
$\ell$-automatic if there is a DFAO
$(\mathcal{S},\\{0,\dots,\ell-1\\},\delta,\\{s_{0}\\},\mathcal{A},\omega)$
such that $u_{n}=\omega(\delta(s_{0},n_{k}\cdots n_{0}))$ for
$(n)_{\ell}=n_{k}\cdots n_{0}$. Similarly, a sequence
$\left(u_{n}\right)_{n\in\mathbb{Z}}$ of elements in $\mathcal{A}$ is
$\ell$-automatic if there is a DFAO
$(\mathcal{S},\\{0,\dots,\ell-1\\},\delta,\\{s_{0},s_{1}\\},\mathcal{A},\omega_{0},\omega_{1})$
such that
$(\mathcal{S},\\{0,\dots,\ell-1\\},\delta,\\{s_{0}\\},\mathcal{A},\omega_{0})$
generates $(u_{n})_{n\geq 0}$ starting at $s_{0}$, and such that if $n<0$,
then $u_{n}=\omega_{1}(\delta(s_{1},\ell-1\,n_{k}\cdots n_{0}))$ for
$(n)_{\ell}=\overline{\ell-1}\,n_{k}\cdots n_{0}$.

If we feed the DFAO starting with the least significant digit, we call this a
_reverse-reading automaton_ ; starting with the most significant non-zero
digit, we call this a _direct-reading automaton_.

Note that in general a DFAO that reads expansions in reverse reading will
generate a different sequence to the same DFAO which reads expansions in
direct reading. Also, an automatic sequence generated by a DFAO in direct
reading can be generated by a (possibly different) DFAO in reverse reading,
and vice versa. The details of the transition from a direct-reading to a
reading-reverse automaton involve reversing all edges in the given automaton;
the resulting automaton may be non-deterministic, i.e., there could be
multiple edges labelled $i$ emanating from one state, and in this case the
automaton needs to be determinised. This is a classical construction, but, it
can blow up the size of the automaton exponentially in $\ell$; see for example
the difference in state size in Figures 1 and 2 which are the direct and
reverse reading automata that generated the (same) substitution in those
examples.

Let $\theta:\,\mathcal{A}^{*}\to\mathcal{A}^{+}$ be a morphism, also called a
_substitution_. Note that the image of any letter is a non-empty word. We will
abuse the notation and write $\theta:\,\mathcal{A}\to\mathcal{A}^{+}$. Let
$\mathbb{I}=\mathbb{N}_{0}$ or $\mathbb{Z}$, according to whether we are to
construct a one-sided or two-sided sequence. Using concatenation, we extend
$\theta$ to $\mathcal{A}^{\mathbb{I}}$. The finiteness of $\mathcal{A}$
guarantees that $\theta$-periodic points, i.e., points x such that
$\theta^{k}(\textit{{x}})=\textit{{x}}$ for some positive $k$, exist. A
substitution is _primitive_ if there exists an integer $n\geq 1$ such that for
all $a,b\in\mathcal{A}$, a occurs in $\theta^{n}(b)$. A substitution $\theta$
is aperiodic if its $\theta$-periodic points are not shift-periodic. In
applications to dynamics we will assume that $\theta$ is primitive and
aperiodic, but our automaton results do not require this constraint.

The following theorem by Cobham [2] connects automatic sequences to
substitution fixed sequences.

######  Theorem 3 (Cobham’s theorem).

A sequence $\textit{{u}}=(u_{n})_{n\geq 0}$ is $\ell$-automatic in direct
reading using $(\mathcal{S},\Sigma_{\ell},\delta,s_{0},\mathcal{A},\omega)$,
where $\delta(s_{0},0)=s_{0}$, if and only if u is the image under a coding of
a right-infinite fixed point of a substitution of constant length $\ell$.

Let $\theta$ be a length-$\ell$ substitution. We write
$\theta(a)=\theta_{0}(a)\cdots\theta_{\ell-1}(a)$; with this notation we see
that for each $0\leq i\leq\ell-1$, we have a map
$\theta_{i}:\mathcal{A}\rightarrow\mathcal{A}$ where, for each
$a\in\mathcal{A}$, $\theta_{i}(a)$ is the $i$-th letter of the word
$\theta(a)$. The map $\theta_{i}$ can be visualised as the $i$-th column of
letters if we stack the substitution words in an
$|\mathcal{A}|\times\ell$-array.

We will assume that our substitution is in simplified form, i.e. $\theta_{0}$
and $\theta_{\ell-1}$ are each idempotents in $\mathcal{A}^{\mathcal{A}}$. One
reason for this is that the automaton will generate the same element $u_{n}$
regardless of the chosen base-$\ell$ expansion of $n$, i.e. the same entry is
obtained if we allow leading zeroes in $(n)_{\ell}$, for $n\in\mathbb{N}_{0}$,
or leading $\ell-1$’s for $n\in-\mathbb{N}$. But also, we can then lightly
strengthen the statement of Cobham’s theorem, so that it concerns bi-infinite
fixed points of substitutions; this is desired as we work with two-sided and
hence invertible shifts. From a dynamical perspective, as we work with
primitive aperiodic substitutions, the dynamical system generated by $\theta$
equals that generated by any of its powers, so there is no loss of generality
in assuming that $\theta$ is in simplified form.

We can then extend the definition of $\theta_{i}$ for $0\leq i\leq\ell-1$ for
a simplified substitution to define $\theta_{n}$ where $n\in\mathbb{Z}$ as
follows. If $n\in\mathbb{N}$ has canonical form $n_{j}\cdots n_{0}$, we define
$\theta_{n}:=\theta_{n_{0}}\circ\dots\circ\theta_{n_{j}}$. Similarly, if
$n\in-\mathbb{N}\backslash\\{-1\\}$ has canonical form $\ell-1\,n_{j}\cdots
n_{0}$, we define
$\theta_{n}:=\theta_{n_{0}}\circ\dots\circ\theta_{n_{j}}\circ\theta_{\ell-1}$.
Finally define $\theta_{-1}:=\theta_{\ell-1}$. We can also extend the notion
of a column map applied to a sequence denoted by (emboldened)
$\boldsymbol{\theta_{r}}:\mathcal{A}^{\mathbb{I}}\to\mathcal{A}^{\mathbb{I}}$,
by setting

(1)
$\boldsymbol{\theta_{r}}\left((u_{n}\right)_{n\in\mathbb{I}})=\left(\theta_{r}(u_{n})\right)_{n\in\mathbb{I}}$

Given a substitution in simplified form, there is a very natural direct
reading DFA that will generate a one-sided fixed point u, given by the proof
of Cobham’s theorem [2]. Its initial state is the first letter $u_{0}$ of u,
its states are labelled by letters in $\mathcal{A}$, and its transition
function is given by the maps $\theta_{i}$, as $\delta(a,i):=\theta_{i}(a)$.
Moreover, although the classical literature concerning automata deals with
one-sided sequences, it is straightforward to modify the statements to
generate two-sided fixed points. Namely, to generate the entry $u_{-n}$ of the
left-infinite fixed point $\dots,u_{-2},u_{-1}$, if
$(n)_{\ell}=\overline{\ell-1}\,n_{j}\cdots n_{0}$, we input the finite word
$\ell-1\,n_{j}\cdots n_{0}$ in direct reading, but starting at the initial
state $u_{-1}$. We summarise some of the information that we will need in the
following lemma; its proof can be found in [2] for the case when
$n\in\mathbb{N}_{0}$, and a simple generalisation works for $n\in-\mathbb{N}$.

######  Lemma 4.

Let $\theta$ be a length-$\ell$ substitution, with fixed point given by the
seed $b\cdot a$. If $n\in\mathbb{N}_{0}$ with $(n)_{\ell}=n_{k}\cdots n_{0}$,
then $u_{n}=\theta_{n_{0}}\circ\cdots\circ\theta_{n_{k}}(a)$. If
$n\in-\mathbb{N}$ with $(n)_{\ell}=(\overline{\ell-1}\,n_{k}\cdots n_{0})$,
then $u_{n}=\theta_{n_{0}}\circ\cdots\circ\theta_{n_{k}}\theta_{\ell-1}(b)$.

######  Example 5.

Define $\theta:\\{a,b,c\\}\rightarrow\\{a,b,c\\}^{3}$ as $a\mapsto acb$,
$b\mapsto baa$, and $c\mapsto bba$. $\theta$ is in its simplified form. A
minimal direct-reading automaton of $\theta$ that generates the fixed point
given by the seed $b\cdot a$ is given in Figure 1. We feed non-negative
integers to the initial state $a$, and negative integers to the initial state
$b$.

$a$$\mathbb{N}_{0}$$b$$-\mathbb{N}$$c$01201,20,12 Figure 1. A minimal direct-
reading automaton for $\theta$ in Example 5

###  2.1. The $\ell$-kernel and the reverse-reading automaton

The _$\ell$ -kernel_ of a sequence $\left(u_{n}\right)_{n\geq 0}$ is the
collection of sequences

$\ker_{\ell}\bigl{(}\left(u_{n}\right)_{n\geq
0}\bigr{)}=\bigl{\\{}\bigl{(}u_{\ell^{e}\cdot n+j}\bigr{)}_{n\geq 0}\,:\,e\geq
0,\,0\leq j\leq\ell^{e}-1\bigr{\\}}.$

The following theorem is due to Eilenberg [1].

######  Theorem 6 (Eilenberg’s theorem).

A sequence u is $\ell$-automatic in reverse reading if and only if it has a
finite $\ell$-kernel. Furthermore, the size of a minimal finite-state
automaton which generates u is the cardinality of $\ker_{\ell}(\textit{{u}})$.

As with Cobham’s theorem, Eilenberg’s theorem is stated for one-sided
sequences. However we can equally define the $\ell$-kernel of a two-sided
sequence, and give a statement and proof of Eilenberg’s theorem that works for
two-sided sequences.

Let $\mathbb{I}=\mathbb{N}_{0}$ or $\mathbb{I}=\mathbb{Z}$. To study the
$\ell$-kernel of a sequence u, we introduce the operators
$\Lambda_{i}:\mathcal{A}^{\mathbb{I}}\rightarrow\mathcal{A}^{\mathbb{I}}$, for
$0\leq i\leq\ell-1$ defined by

$\Lambda_{i}(\left(u_{n}\right)_{n\in\mathbb{I}}):=\left(u_{\ell\cdot
n+i}\right)_{n\in\mathbb{I}}.$

In the case where the sequence u takes values in a finite field and
$\ell=p^{j}$, with $p$ prime, and when $\mathbb{I}=\mathbb{N}$, the operators
$\Lambda_{i}$ are called the Cartier operaters, except that they act on formal
power series instead of sequences.

The usual construction of a reverse-reading automaton that generates a one-
sided automatic sequence, e.g., the one described in [1], has states labelled
with elements of $\ker_{\ell}\bigl{(}\textit{{u}}\bigr{)}$, and output
obtained by projecting the 0-indexed entry, $\pi_{0}((u_{n})):=u_{0}$. For v,
a sequence labelling one of the states, the transition function $\delta$ is
$\delta(\textit{{v}},i):=\Lambda_{i}(\textit{{v}})$. We will give an alternate
construction of the reverse-reading automaton which will be useful for our
application to Toeplitz sequences. To do this, we link the operators
$\Lambda_{i}$ to the column maps ${\bf\theta}_{i}$. Recall the definition of
$\boldsymbol{\theta_{r}}$ in (1).

######  Lemma 7.

Let $\textit{{u}}=\left(u_{n}\right)_{n\in\mathbb{I}}$ be a fixed point of a
length-$\ell$ substitution $\theta$. Then for each $0\leq r\leq\ell-1$, we
have

(2) $\Lambda_{r}(\textit{{u}})=\boldsymbol{\theta_{r}}(\textit{{u}}).$

###### Proof.

Given a fixed point u, note that by definition, for each $k\in\mathbb{I}$, we
have $u_{k\cdot\ell},\dots,u_{k\cdot\ell+\ell-1}=\theta(u_{k})$. we have
$\theta(\textit{{u}})=\textit{{u}}$ by definition. Therefore
$u_{k\cdot\ell+r}=\theta_{r}(u_{k})$, so that
$\boldsymbol{\theta_{r}}(\textit{{u}})=(\theta_{r}(u_{k}))_{k\in\mathbb{I}}=\left(u_{k\cdot\ell+r}\right)_{k\in\mathbb{I}}$.
On the other hand,

$\Lambda_{r}(\textit{{u}})=\Lambda_{r}(\left(u_{k}\right)_{k\in\mathbb{I}})=\left(u_{k\cdot\ell+r}\right)_{k\in\mathbb{I}},$

from which the result follows. ∎

As with the definition of the column maps $\theta_{n}$, we can extend the
defintion of $\Lambda_{i}$ for $0\leq i\leq\ell-1$ to the definition of
$\Lambda_{n}$, where $n\in\mathbb{N}$ as follows: if $(n)_{\ell}=n_{j}\cdots
n_{0}$, with $n_{j}\neq 0$, then
$\Lambda_{n}:=\Lambda_{n_{j}}\circ\dots\circ\Lambda_{n_{0}}$. Similarly, if
$n\in-\mathbb{N}$ and $(n)_{\ell}=\overline{\ell-1}\,n_{j}\cdots n_{0}$ with
$n_{j}\neq\ell-1$, then
$\Lambda_{n}:=\Lambda_{\ell-1}\circ\Lambda_{n_{j}}\circ\dots\circ\Lambda_{n_{0}}$.
We have the following as a corollary of Lemma 7.

######  Corollary 8.

Let $\textit{{u}}=\left(u_{n}\right)_{n\in\mathbb{I}}$ be a fixed point of a
length-$\ell$ substitution $\theta$. Then for any $n_{0},n_{1},\dots,n_{d}$
with $0\leq n_{i}\leq\ell-1$ for each $i$, we have

(3)
$\Lambda_{n_{d}}\circ\Lambda_{n_{d-1}}\circ\dots\circ\Lambda_{n_{0}}(\textit{{u}})=\boldsymbol{\theta_{n_{0}}}\circ\boldsymbol{\theta_{n_{1}}}\circ\dots\circ\boldsymbol{\theta_{n_{d}}}(\textit{{u}})$

###### Proof.

Given $d\in\mathbb{N}$ and $0\leq n_{0},n_{1},\dots,n_{d}\leq\ell-1$, it can
be shown as in the proof of Lemma 7 by induction that

$({\boldsymbol{\theta}}_{n_{0}}\circ{\boldsymbol{\theta}}_{n_{1}}\circ\dots\circ{\boldsymbol{\theta}}_{n_{d}}(\textit{{u}}))_{n}=u_{n\cdot\ell^{d+1}+n_{d}\cdot\ell^{d}+\dots+n_{1}\cdot\ell+n_{0}}$

for each $n\in\mathbb{I}$. We shall show that

(4)
$\left(\Lambda_{n_{d}}\circ\dots\circ\Lambda_{n_{1}}\circ\Lambda_{n_{0}}({u_{n}})\right)_{n\in\mathbb{I}}=\left(u_{n\cdot\ell^{d+1}+n_{d}\cdot\ell^{d}+\dots+n_{1}\cdot\ell+n_{0}}\right)_{n\in\mathbb{I}}$

which will prove the corollary.

If $d=0$, the assertion follows by Lemma 7. Inductively, assume that
$\left(\Lambda_{n_{d-1}}\circ\dots\circ\Lambda_{n_{1}}\circ\Lambda_{n_{0}}({u_{n}})\right)_{n\in\mathbb{I}}=\left(u_{n\cdot\ell^{d}+n_{d-1}\cdot\ell^{d-1}+\dots+n_{1}\cdot\ell+n_{0}}\right)_{n\in\mathbb{I}}$
holds for $d\in\mathbb{N}$. Then

$\displaystyle(\Lambda_{n_{d}}\circ\dots\circ\Lambda_{n_{1}}\circ\Lambda_{n_{0}}({u_{n}}))_{n\in\mathbb{I}}$
$\displaystyle=(\Lambda_{n_{d}}(\Lambda_{n_{d-1}}\circ\dots\circ\Lambda_{n_{1}}\circ\Lambda_{n_{0}}({u_{n}})))_{n\in\mathbb{I}}$
$\displaystyle=\Lambda_{n_{d}}\left(u_{n\cdot\ell^{d}+n_{d-1}\cdot\ell^{d-1}+\dots+n_{1}\cdot\ell+n_{0}}\right)_{n\in\mathbb{I}}$
$\displaystyle=\left(u_{n\cdot{\ell^{d+1}}+n_{d}\cdot\ell^{d}+{n_{d-1}}\cdot{\ell^{d-1}}+\dots+{n_{1}}\cdot\ell+{n_{0}}}\right)_{n\in\mathbb{I}}$

and the induction is complete. ∎

Let $\operatorname{id}:\mathcal{A}\rightarrow\mathcal{A}$ denote the identity
map.

######  Definition 9.

[The semigroup of a substitution] Let $\theta$ be a length-$\ell$ substitution
on $\mathcal{A}$ in simplified form. The semigroup of $\theta$, denoted
$\mathcal{S}_{\theta}$ is the semigroup in $\mathcal{A}^{\mathcal{A}}$ defined
by

$\mathcal{S}_{\theta}:=\bigcap_{n\geq
1}\langle\operatorname{id},{\theta^{n}}_{i}:0\leq i\leq\ell^{n}-1\rangle.$

Note that $\langle\operatorname{id},{\theta^{n}}_{i}:0\leq
i\leq\ell^{n}-1\rangle\subset\langle\operatorname{id},{\theta^{d}}_{i}:0\leq
i\leq\ell^{d}-1\rangle$ whenever $d|n$, so there is an $n$ such that
$\mathcal{S}_{\theta}=\langle\operatorname{id},{\theta^{n}}_{i}:0\leq
i\leq\ell^{n}-1\rangle$. Thus by working with power of $\theta$, if necessary,
we can assume that
$\mathcal{S}_{\theta}:=\langle\operatorname{id},{\theta}_{i}:0\leq
i\leq\ell-1\rangle.$

######  Example 10.

Define $\theta:\\{a,b\\}\rightarrow\\{a,b\\}^{2}$ as $a\mapsto ab$ and
$b\mapsto aa$. This is the period-doubling substitution, and its simplified
form is $\theta^{2}$; so $a\mapsto abaa$ and $b\mapsto abab$. we have for
$n=2$,

$\langle\operatorname{id},{\theta^{2}}_{i}:0\leq i\leq
3\rangle=\\{\operatorname{id},(a,a)^{T},(b,b)^{T}\\}\subsetneq\\{\operatorname{id},(a,a)^{T},(b,b)^{T},(b,a)^{T}\\}=\langle\operatorname{id},{\theta}_{i}:0\leq
i\leq 1\rangle.$

It can be verified that $\langle\operatorname{id},{\theta^{2}}_{i}:0\leq i\leq
3\rangle=\langle\operatorname{id},{\theta^{n}}_{i}:0\leq
i\leq\ell^{n}-1\rangle$ for each $n\geq 2$, so that

$\mathcal{S}_{\theta}=\bigcap_{n\geq
1}\langle\operatorname{id},{\theta^{n}}_{i}:0\leq
i\leq\ell^{n}-1\rangle=\\{\operatorname{id},(a,a)^{T},(b,b)^{T}\\}.$

If $\theta$ is bijective, i.e., if each $\theta_{i}$ is a bijection, then
$\mathcal{S}_{\theta}$ is a subgroup. Its importance is recognized in [8], and
it is related to the Ellis semigroup in [7]. In general the semigroup
generated by the maps $\\{\theta_{0},\dots,\theta_{\ell-1}\\}$ does not
necessarily contain $\operatorname{id}$.

The proof of the following lemma is a straightforward induction argument.

######  Lemma 11.

If $\pi_{0}:\mathcal{A}^{\mathbb{I}}\rightarrow\mathcal{A}$ is the map which
projects a sequence $(u_{n})$ to its zero-indexed entry $u_{0}$, and if
$n\in\mathbb{N}_{0}$ with $(n)_{\ell}=(n_{k}\cdots n_{0})$ its canonical
representation, then
$u_{n}=\pi_{0}\left(\Lambda_{n_{k}}\circ\cdots\circ\Lambda_{n_{0}}({\textit{{u}}})\right)$.
Similarly, if $\pi_{-1}:\mathcal{A}^{\mathbb{Z}}\rightarrow\mathcal{A}$ is
defined as $\pi_{-1}((u_{n})_{n\in Z}):=u_{-1}$, and $n\in-\mathbb{N}$ with
$(n)_{\ell}=(\overline{\ell-1}\,n_{k}\cdots n_{0})$, then
$u_{n}=\pi_{-1}\left(\Lambda_{\ell-1}\circ\Lambda_{n_{k}}\circ\cdots\circ\Lambda_{n_{0}}({\textit{{u}}})\right)$.

We write elements $s\in\mathcal{A}^{\mathcal{A}}$ as vectors, i.e., if
$\mathcal{A}:=\\{a_{0},\dots,a_{d}\\}$, then
$s=\begin{pmatrix}b_{0},\dots,b_{d}\end{pmatrix}^{T}$ if and only if
$s(a_{i})=b_{i}$ for each $i$.

######  Theorem 12.

Let $\theta$ be a length-$\ell$ substitution on $\mathcal{A}$, in simplified
form. Then there is a minimal automaton
$\mathcal{M}=(\mathcal{S},\Sigma,\delta,s_{0},s_{1})$, whose set of states
$\mathcal{S}=\mathcal{S}_{\theta}$, the structure semigroup of $\theta$, such
that for any fixed point of $\theta$ with seed ${\textit{{u}}}=u_{l}\cdot
u_{r}$, there are output maps
$\omega_{l},\omega_{r}:\mathcal{S}\rightarrow\mathcal{A}$ such that
$\mathcal{M},\omega_{l},\omega_{r}$ generates u in reverse reading.

###### Proof.

Suppose that $\theta$ is defined on $\\{a_{0},\dots,a_{j}\\}$. We will assume
that we are given a bi-infinite fixed point, generated by $a_{l}\cdot a_{r}$,
where $\theta(a_{l})$ ends with $a_{l}$ and $\theta(a_{r})$ starts with
$a_{r}$.

We first define the finite-state automaton
$\mathcal{M}=(\mathcal{S},\Sigma,\delta,s_{0},\mathcal{A},\omega_{r})$ that
generates the right part of u. We let $\mathcal{S}=\mathcal{S}_{\theta}$, and
we will abuse notation and henceforth think of states as functions, even if
technically they are labelled by functions. We let
$\Sigma=\\{0,\dots,\ell-1\\}$, $\mathcal{A}:=\\{a_{0},\dots,a_{j}\\}$ and we
let the initial state $s_{0}:=\operatorname{id}$. The output map
$\omega_{r}:\mathcal{S}_{\theta}\rightarrow\mathcal{A}$ will be
$p_{a_{r}}(s)$, which denotes projection to $s(a_{r})$. It remains to define
the transition map $\delta$.

For $0\leq i\leq\ell-1$, define $\delta(\operatorname{id},i):=\theta_{i}$.
Note that by definition any element $s$ in $\mathcal{S}_{\theta}$ can be
written as $s=\theta_{n_{0}}\circ\dots\circ\theta_{n_{k}}$ for some $0\leq
n_{0},\dots,n_{k}\leq\ell-1$. Now set
$\delta(s,i):=\theta_{n_{0}}\circ\cdots\circ\theta_{n_{k}}\theta_{i}(s_{0})$;
this is well-defined.

It remains to show that $\mathcal{M}$ generates the right-infinite part of u.
If $(n)_{\ell}=n_{k}\cdots n_{0}$, then

$u_{n}\stackrel{{\scriptstyle
Lemma\,\,\ref{lem:projection}}}{{=}}\pi_{0}\left(\Lambda_{n_{k}}\circ\cdots\circ\Lambda_{n_{0}}({\textit{{u}}})\right)\stackrel{{\scriptstyle
Corollary\,\,\ref{cor:column=cartier}}}{{=}}\pi_{0}\left({\boldsymbol{\theta}}_{n_{0}}\circ\cdots\circ{\boldsymbol{\theta}}_{n_{k}}({\boldsymbol{u}})\right)\stackrel{{\scriptstyle
Lemma\,\,\ref{lem:projection-0}}}{{=}}{\theta}_{n_{0}}\circ\cdots\circ{\theta}_{n_{k}}(p_{a_{r}}(\operatorname{id})),$

and
${\theta}_{n_{0}}\circ\cdots\circ{\theta}_{n_{k}}(p_{a_{r}}(\operatorname{id}))=p_{a_{r}}\left({\theta}_{n_{0}}\circ\cdots\circ{\theta}_{n_{k}}(\operatorname{id})\right)$
as $p_{a_{r}}$ commutes with each $\theta_{i}$. The result follows.

The left part of ${\boldsymbol{u}}$ will be generated similarly, except with
the output map $\omega_{a_{l}}$. Furthermore, the states are in one-to-one
correspondence with the kernel of $\theta$. In particular, if a state is
labelled by $s\in\mathcal{S}_{\theta}$, then it represents that sequence ${\bf
s}({\textit{{u}}})$. Thus, $\mathcal{M}$ is minimal by Eilenberg’s theorem.

∎

######  Example 13.

We look at Example 5 again. Note that since for each $n$,
${\theta^{n}}_{0}=(a,b,b)^{T}$, ${\theta^{n}}_{2}=(b,a,a)^{T}$, and
${\theta^{n}}_{\alpha}\in\langle(c,a,b)\rangle$ for
$\alpha=\sum_{i=0}^{n-1}3^{i}$, we have that
$\langle\operatorname{id},{\theta^{n}}_{i},0\leq
i\leq\ell^{n}-1\rangle=\langle\operatorname{id},{\theta}_{i},0\leq i\leq
2\rangle$, so
$\mathcal{S}_{\theta}=\langle\operatorname{id},{\theta}_{i},0\leq i\leq
2\rangle$. A minimal reverse-reading automaton for $\theta$, with states
labelled using $\mathcal{S}_{\theta}$, is given in Figure 2.

Figure 2. Reverse-reading automaton for Example 13; any of the intersecting
edges with the corresponding labels are coloured differently for more clarity.

##  3\. Toeplitz sequences

There are many definitions of Toeplitz sequences, and a vast literature on
Toeplitz shifts; we use the notions which bring us directly to our setting,
which is that of substitution sequences.

###  3.1. The column number of a substitution

Let $\theta$ be a primitive length-$\ell$ substitution with fixed point u, and
with $(X_{\theta},\sigma)$ infinite. The height $h=h(\theta)$ of $\theta$ is
defined as

$h(\theta):=\max\\{n\geq 1:\gcd(n,\ell)=1,n|\gcd\\{a:u_{a}=u_{0}\\}\\}\,.$

If $h>1$, this means that $\mathcal{A}$ decomposes into $h$ disjoint subsets:
$\mathcal{A}_{1}\cup\ldots\cup\mathcal{A}_{h}$, where a symbol from
$\mathcal{A}_{i}$ is always followed by a symbol from
$\mathcal{A}_{i+1\mod{h}}$ [4]. Also, such a system is a constant height
suspension over another length-$\ell$ substitution shift which has trivial
height.

######  Theorem 14.

[4, Remark 9, Lemmas 17 and 19] If $\theta$ is a primitive aperiodic
length-$\ell$ substitution defined on $\mathcal{A}$, and has height $h$, then
there is a length-$\ell$ substitution $\bar{\theta}$ which has trivial height,
and $(X_{\theta},\sigma)\cong(X_{\bar{\theta}}\times\\{0,\ldots,h-1\\},T)$
where

$T(x,i)\coloneqq\begin{cases}(x,i+1)&\text{ if }0\leq i<h-1\\\
(\bar{\sigma}(x),0)&\text{ if }i=h-1\end{cases}$

The substitution $\bar{\theta}$ is called the pure base of $\theta$. If
$\theta$ has height $h=1$, then it is its own pure base.

Let $\theta$ have pure base $\bar{\theta}$. We say that $\theta$ has column
number $c$. if for some $k\in\mathbb{N}$, and some $(i_{1},\ldots,i_{k})$,
$|\bar{\theta}_{i_{1}}\ldots\bar{\theta}_{i_{k}}(\mathcal{A})|=c$ and $c$ is
the least such number. In particular, if $\theta$ has column number one, then
we will say that $\theta$ has a coincidence.

###  3.2. Odometers

Given a sequence $(\ell_{n})$ of natural numbers, we work with the group

$\mathbb{Z}_{(\ell_{n})}:=\prod_{n}\mathbb{Z}/\ell_{n}\mathbb{Z},$

where the group operation is given by coordinate-wise addition with carry. For
a detailed exposition of equivalent definitions of $\mathbb{Z}_{(\ell_{n})}$,
we refer the reader to [5]. Endowed with the product topology over the
discrete topology on each $\mathbb{Z}/\ell_{n}\mathbb{Z}$, the group
$\mathbb{Z}_{(\ell_{n})}$ is a compact metrizable topological group, where the
unit $z=\dots 001$, which we simply write as $z=1$, is a topological
generator. We write elements $(z_{n})$ of $\mathbb{Z}_{(\ell_{n})}$ as left-
infinite sequences $\dots z_{2}z_{1}$ where
$z_{n}\in\mathbb{Z}/\ell_{n}\mathbb{Z}$, so that addition in
$\mathbb{Z}_{(\ell_{n})}$ has the carries propagating to the left as usual in
$\mathbb{Z}$. If $\ell_{n}=\ell$ is constant, then
$\mathbb{Z}_{(\ell_{n})}=\mathbb{Z}_{\ell}$ is the classical ring of
$\ell$-adic integers.

With the above notation, an odometer is a dynamical system
$(\mathbb{Z}_{(\ell_{n})},+1)$. Odometers are _equicontinuous_ , i.e., the
family $\\{+n:n\in\mathbb{Z}\\}$ is equicontinuous.

An equicontinuous factor of $(X,T)$ is maximal if any other equicontinuous
factor of $(X,T)$ factors through it; the maximal equicontinuous factor always
exists. A Toeplitz shift is a symbolic shift $(X,\sigma)$,
$X\subset\mathcal{A}^{\mathbb{Z}}$, which is a somewhere one-to-one extension
of an odometer. The following theorem tells us that the Toeplitz length-$\ell$
substitution shifts are precisely those that have column number one.

######  Theorem 15.

Let $\theta$ be an aperiodic primitive, length-$\ell$ substitution of height
$h$. Then the maximal equicontinuous factor of $(X_{\theta},\sigma)$ is the
odometer $(\mathbb{Z}_{\ell}\times\mathbb{Z}/h\mathbb{Z},(+1,+1))$.
Furthermore $(X_{\theta},\sigma)$ is a Toeplitz shift if and only if $\theta$
has column number one.

###  3.3. Toeplitz sequences

A Toeplitz sequence is a sequence that is obtained by filling its entries one
arithmetic sequence of indices at a time, as follows. Let
$\mathbb{I}=\mathbb{N}_{0}$ or $\mathbb{Z}$. Let $\mathcal{A}$ be a finite
alphabet. For $\textit{{x}}\in\mathcal{A}^{\mathbb{I}}$, $k\in\mathbb{N}$ and
$a\in\mathcal{A}$, we define

$\mbox{Per}_{k}(\textit{{x}},a):=\\{n\in\mathbb{I}:x_{n^{\prime}}=a\mbox{ for
all }n^{\prime}\equiv n\bmod k\\},$
$\mbox{Per}_{k}(\textit{{x}}):=\bigcup_{a\in\mathcal{A}}\mbox{Per}_{k}(\textit{{x}},a),$

and

$\mbox{Aper}(\textit{{x}}):=\mathbb{I}\backslash\left(\bigcup_{k\in\mathbb{N}}\mbox{Per}_{k}(\textit{{x}})\right).$

We say that the sequence x is Toeplitz if
$\mbox{Aper}(\textit{{x}})=\emptyset$. It can be shown that a minimal shift
$(X,\sigma)$ is Toeplitz if and only if contains a Toeplitz sequence. Note
that a substitution shift may be Toeplitz, but its fixed points are not
necessarily Toeplitz sequences, as the following example shows.

######  Example 16.

Recall the period-doubling substitution $\theta$ from Example 10. A minimal
reverse-reading automaton that generates the period-doubling sequence x with
states as $\theta_{i}\in\mathcal{A}^{\mathcal{A}}$ is the following:

$\begin{matrix}a\\\ b\end{matrix}$$\begin{matrix}a\\\
a\end{matrix}$$\begin{matrix}b\\\ b\end{matrix}$0,2130,1,2,30,1,2,3 Figure 3.
A minimal reverse-reading automaton generating x

Now, we generate the fixed point of $\theta^{2}$ using the algorithm used to
define the Toeplitz sequence. To generate the sequence
$\textit{{x}}\in\mathcal{A}^{\mathbb{Z}}$, we take $k=2^{i}$ for each
iteration of choice of an arithmetic sequence; i.e., $i=1,2,\dots$ and fill
the letters $a$ and $b$ on alternate iterations starting with $a$. For $i=1$,
$\mbox{Per}_{2}(\textit{{x}},a)=\\{2\mathbb{I}\\}$;
$i=2\implies\mbox{Per}_{4}(\textit{{x}},b)=\\{4\mathbb{I}+1\\}$;
$i=3\implies\mbox{Per}_{8}(\textit{{x}},a)=\\{8\mathbb{I}+3\\}$;
$i=4\implies\mbox{Per}_{16}(\textit{{x}},b)=\\{16\mathbb{I}+7\\}$, and so on.
Then $\mbox{Per}_{k}(\textit{{x}})=\\{\mathbb{I}\\}\backslash\\{-1\\}$. So,
$\mbox{Aper}(\textit{{x}})=\\{-1\\}\neq\emptyset$. Another way to realise this
is that every entry in the two fixed points possible is the same except
$x_{-1}$. Hence, the fixed points are not Toeplitz sequences. The column
number of the substitution $\theta$ is one, which confirms that it generates a
Toeplitz shift. Although, notice that the right-infinite fixed point of
$\theta$ is in fact a Toeplitz sequence.

Recall that each state $s$ in the semigroup automaton $\mathcal{M}$ is
labelled by a map $f_{s}:\mathcal{A}\rightarrow\mathcal{A}$. Call the state
$s$ a $k$-vertex if $|k|=|\mbox{Im}f_{s}|$. If $s$ is a $1$-vertex, then
$f_{s}=\pi_{a}$, the projection to the letter $a$ for some $a$.

######  Definition 17 (Reduced graph for a substitution).

Let $\theta$ have a coincidence. If we remove from its semigroup automaton all
$1$-vertices, and all edges leading to them, we are left with the reduced
graph for $\theta$.

From Example 16, the reduced graph for $\theta^{2}$ would be the following:

$\begin{matrix}a\\\ b\end{matrix}$3 Figure 4. The reduced graph for
$\theta^{2}$

This confirms again, that the aperiodic part of the sequence is $-1$, as the
infinite loop $\stackrel{{\scriptstyle\xleftarrow{\quad\qquad}}}{{\dots 333}}$
is the base-$4$ representation of $-1$. If $(n)_{\ell}=n_{j-1}\dots n_{0}$, we
write $|n|_{\ell}=j$.

######  Theorem 18.

Let $\theta$ be a primitive aperiodic length-$\ell$ substitution in simplified
form, of trivial height and let $\mathcal{M}$ be the semigroup finite-state
automaton which generates $\theta$ in reverse reading. Let u be an bi-infinite
fixed point of $\theta$. If $(X_{\theta},\sigma)$ is Toeplitz, then

$\mbox{Per}(\textit{{u}})=\\{n\in\mathbb{Z}:\delta(u_{0},(n)_{\ell})=\pi_{u_{n}}\mbox{
and }\delta(u_{-1},(n-p^{|n|_{\ell}})_{\ell})=\pi_{u_{n}}\\}.$

###### Proof.

First suppose that $n\in\mathbb{N}_{0}$ and that
$n\in\mbox{Per}(\textit{{u}})$. As $n\in\mbox{Per}(\textit{{u}})$, we have
$(u_{n+kp^{j}})_{k\in\mathbb{Z}}$ is constant sequence, i.e.
$\Lambda_{n_{j-1}}\dots\Lambda_{n_{0}}(\textit{{u}})$ is constantly equal to
$a=u_{n}$, so that $(u_{n+k\cdot p^{j}})_{k\in\mathbb{N}}$ is the constant
sequence $\bar{a}$. By the construction of $\mathcal{M}$, we have
$\delta(u_{0},(n)_{\ell})=\pi_{a}$, the projection to the letter $a$.
Similarly, $(u_{n+k\cdot p^{j}})_{k\in-\mathbb{N}}=\bar{a}$, so in particular
letting $k=-1$, $\delta(u_{-1},(n-p^{j})_{\ell})=\pi_{a}$.

Conversely, suppose that $n$ is such that $u_{n}=a$, and also
$\delta(u_{0},(n)_{\ell})=\pi_{a}$ and
$\delta(u_{-1},(n-p^{|n|_{\ell}})_{\ell})=\pi_{a}$. Let $j=|n|_{\ell}$. First
suppose that $k\in\mathbb{N}_{0}$. Then

$\delta(u_{0},(n+k\cdot
p^{j})_{\ell})=\theta_{n}\circ\theta_{k}\left(\begin{pmatrix}a_{1},a_{2},\dots,a_{d}\end{pmatrix}^{T}\right)=\pi_{a}\circ\theta_{k}=\pi_{a}.$

Similarly, if $k\in-\mathbb{N}$, then
$\delta(u_{-1},(n-p^{j}+(k+1)p^{k})_{\ell})=\theta_{n-p^{j}}\circ\theta_{k+1}\left(\begin{pmatrix}a_{1},a_{2},\dots,a_{d}\end{pmatrix}\right)=\pi_{a}\circ\theta_{n-p^{j}}=\pi_{a}$.

∎

## References

  * [1] Jean-Paul Allouche and Jeffrey Shallit, _Automatic sequences_ , Cambridge University Press, Cambridge, 2003, Theory, applications, generalizations. MR 1997038
  * [2] Alan Cobham, _Uniform tag sequences_ , Math. Systems Theory 6 (1972), 164–192. MR 457011
  * [3] Ethan M. Coven, Anthony Quas, and Reem Yassawi, _Computing automorphism groups of shifts using atypical equivalence classes_ , Discrete Anal. (2016), Paper No. 3, 28. MR 3533302
  * [4] F. M. Dekking, _The spectrum of dynamical systems arising from substitutions of constant length_ , Z. Wahrscheinlichkeitstheorie und Verw. Gebiete 41 (1977/78), no. 3, 221–239. MR 461470
  * [5] Tomasz Downarowicz, _Survey of odometers and Toeplitz flows_ , Algebraic and topological dynamics, Contemp. Math., vol. 385, Amer. Math. Soc., Providence, RI, 2005, pp. 7–37. MR 2180227
  * [6] Konrad Jacobs and Michael Keane, _$0-1$ -sequences of Toeplitz type_, Z. Wahrscheinlichkeitstheorie und Verw. Gebiete 13 (1969), 123–131. MR 255766
  * [7] Johannes Kellendonk and Reem Yassawi, _The Ellis semigroup of bijective substitutions_ , Groups Geom. Dyn. 16 (2022), no. 1, 29–73. MR 4424962
  * [8] Mariusz Lemańczyk and Mieczysław K. Mentzen, _On metric properties of substitutions_ , Compositio Math. 65 (1988), no. 3, 241–263. MR 932072
  * [9] John C. Oxtoby, _Ergodic sets_ , Bull. Amer. Math. Soc. 58 (1952), 116–136. MR 47262
  * [10] Susan Williams, _Toeplitz minimal flows which are not uniquely ergodic_ , Z. Wahrsch. Verw. Gebiete 67 (1984), no. 1, 95–107. MR 756807