# Multiqubit state learning with entangling quantum generative adversarial
networks

S. E. Rasmussen<EMAIL_ADDRESS>Department of Physics and Astronomy, Aarhus
University, DK-8000 Aarhus C, Denmark  Department of Chemistry, Aarhus
University, DK-8000 Aarhus C, Denmark  Kvantify Aps, DK-2300 Copenhagen S,
Denmark  N. T. Zinner<EMAIL_ADDRESS>Department of Physics and Astronomy,
Aarhus University, DK-8000 Aarhus C, Denmark  Kvantify Aps, DK-2300 Copenhagen
S, Denmark

###### Abstract

The increasing success of classical generative adversarial networks (GANs) has
inspired several quantum versions of GANs. Fully quantum mechanical
applications of such quantum GANs have been limited to one- and two-qubit
systems. In this paper, we investigate the entangling quantum GAN (EQ-GAN) for
multiqubit learning. We show that the EQ-GAN can learn a circuit more
efficiently compared with a SWAP test. We also consider the EQ-GAN for
learning eigenstates that are variational quantum eigensolver
(VQE)-approximated, and find that it generates excellent overlap matrix
elements when learning VQE states of small molecules. However, this does not
directly translate into a good estimate of the energy due to a lack of phase
estimation. Finally, we consider random state learning with the EQ-GAN for up
to six qubits, using different two-qubit gates, and show that it is capable of
learning completely random quantum states, something which could be useful in
quantum state loading.

Generative machine learning has been highly influenced by generative
adversarial networks (GANs) Goodfellow  _et al._ (2014). In the past couple of
years, GAN applications have increased rapidly; from generating photorealistic
images, Karras  _et al._ (2019) and improving video game resolution Wang  _et
al._ (2018) to modeling dark matter Mustafa  _et al._ (2019) and improving
astrophysical images Schawinski  _et al._ (2017).

The idea behind GANs is that two networks play a game where one network, the
generator, is trying to fool the other network, the discriminator, by
producing realistic data. Thus, the generator is trained indirectly through
the discriminator, hence the name adversarial. Training GANs can be formidably
difficult with vanishing gradients and mode collapse Arjovsky  _et al._ (2017)
and convergence problems Salimans  _et al._ (2016); Mescheder  _et al._
(2018).

The success of classical GANs has inspired quantum physicists to develop a
quantum version, a so-called QuGAN. This was first proposed by the authors of
Refs. Dallaire-Demers and Killoran (2018); Lloyd and Weedbrook (2018), who
argued that a fully quantum mechanical GAN reaches its Nash equilibrium
Osborne and Rubinstein (1994) when the data are reproduced correctly. This
QuGAN does, however, not always converge. In certain cases, it oscillated
between a finite set of states due to mode collapse, as it suffers from a non
unique Nash equilibrium Niu  _et al._ (2022).

In the original proposal, the QuGAN was shown to be able to reproduce a
controlled-not (cnot) gate Dallaire-Demers and Killoran (2018). QuGANs have
been used to approximate pure quantum states Benedetti  _et al._ (2019) and
the modified National Institute of Standards and Technology (MNIST) datasets
Huang  _et al._ (2021), and a proof-of-principle experiment has shown their
feasibility Hu  _et al._ (2019). QuGANs have also been developed with a
quantum generator and a classical discriminator and used on classical data
sets, i.e., financial modeling Zoufal  _et al._ (2019), discrete and
continuous distributions Situ  _et al._ (2020); Romero and Aspuru-Guzik
(2019), the Bars-and-Stripes data set Zeng  _et al._ (2019), and sampling
particle traces Chang  _et al._ (2021). Furthermore, they have been discussed
using a quantum generator and a SWAP test discriminator Stein  _et al._
(2021). In this paper, we consider a fully entangling quantum GAN (EQ-GAN),
i.e., generator, discriminator, and data are quantum mechanical. In its
proposal, the EQ-GAN is shown to be able to generate a single-qubit gate Niu
_et al._ (2022). We apply this EQ-GAN to multiqubit states of up to six
generator qubits and 12 discriminator qubits. We apply it to states found
using a variational quantum eigensolver (VQE) algorithm, to see whether it can
reproduce data from another quantum algorithm and to completely random quantum
states. Generation of quantum states has previously been studied using
classical neural networks Carleo and Troyer (2017); Carleo  _et al._ (2017);
Liu2019; Hibat-Allah  _et al._ (2020); Che  _et al._ (2022).

The paper is structured as follows: Section I presents a short introduction to
the EQ-GAN and the circuit architecture. In Section II we explain how the
simulation is performed, and then in Section II.1 we discuss results for
random two-qubit circuits, in Section II.2 we use an EQ-GAN for learning VQE
states, and in Section II.3 we discuss results when the EQ-GAN is learning
completely random multiqubit states. Finally, we present our conclusion and
the outlook.

##  I EQ-GAN

Figure 1: Parametrized (destructive) SWAP test used as discriminator in this
paper. In order to obtain the original SWAP test one should replace the
controlled-z (cz) gates by a cnot gate, and remove the $R_{z}$ gates and the
Hadamard gates on the qubits representing the $\ket{\psi}$. This is equivalent
to setting $\theta_{i}=0$ for all $i$. Measuring the state of each qubit after
the SWAP test we take a bit wise AND operation of the output $O_{1}$ and
$O_{2}$. If this yields $-1$, then the two states are identical. (a) The one-
qubit case considered in the original proposal Niu  _et al._ (2022). (b) The
multiqubit case used in this paper.

A classical GAN consists of two neural networks: a generator and a
discriminator. The discriminator randomly receives data either from the real
world, denoted as real data, or from the generator, denoted as generated or
fake data. The data can be anything. The discriminator’s task is to determine
whether the data it receives are real or fake. On the other hand, the
generator is trained to generate data to fool the discriminator. The two
networks are thus adversaries: one trying to minimize the cost function and
the other trying to maximize it. For a properly trained GAN, equilibrium is
reached when the discriminator can only randomly guess whether the data are
real or fake. For the details of the cost function see Appendix A.

In the original proposal of a QuGAN Dallaire-Demers and Killoran (2018); Lloyd
and Weedbrook (2018) it was proposed to exchange the neural networks with
parametrized quantum circuits (PQCs) also known as quantum neural networks
Gupta and Zia (2002); Farhi and Neven (2018); Beer  _et al._ (2020). Here, the
discriminator becomes a quantum circuit which is then given a real quantum
state or one generated from another quantum circuit.

This paper considers the entangling quantum GAN (EQ-GAN) Niu  _et al._ (2022).
In the EQ-GAN, real and fake (generated) quantum states are passed to the
discriminator simultaneously such that the discriminator is allowed to
entangle these. The real data could be some quantum state measured in an
experiment, or it could come from some unknown quantum circuit. We denote the
unitary producing this state $R$; it may or may not be dependent on a set of
parameters. The fake quantum state is generated by a PQC given by some circuit
Ansatz. We denote the generator $G(\bm{\theta}_{g})$ where $\bm{\theta}_{d}$
denotes the parameters of the generator.

The discriminator can be either untrainable as a SWAP test Buhrman  _et al._
(2001); Garcia-Escartin and Chamorro-Posada (2013); Stein  _et al._ (2021) or
a trainable circuit capable of learning the SWAP test. The authors of Ref. Niu
_et al._ (2022) show that the trainable discriminator outperforms the perfect
SWAP test in the presence of noise. In this paper, we consider both cases and
denote the discriminator $D(\bm{\theta}_{d})$ where $\bm{\theta}_{d}$ is the
discriminator’s parameters. The discriminator circuit can be seen in Fig. 1(b)
for the multiqubit case considered in this paper. It should be compared with
the single-qubit case in Fig. 1(a) which was considered in the original
proposal Niu  _et al._ (2022).

##  II Results

Figure 2: Example of one layer of the circuit Ansatz,
$G_{l}(\bm{\theta}_{l})$, for four qubits. First, an Euler rotation is
performed on each qubit, followed by a coupling using the given entangling
gate. Here it is shown using fixed cnot gates.

We simulate our EQ-GAN following the approach described in Ref. Broughton  _et
al._ (2021). We prepare our circuits using cirq and then add them as layers
using tensortlow quantum. We compile the model using tensorflow and optimize
it using the adam optimizer, with learning rates of 0.01. We introduce noise
in our system on the rotation angle of all of the gates. We do this by adding
after each gate another gate with a random error. On the single-qubit
rotations we add an error to the rotation angle sampled from a normal
distribution with $\mu=0.06$ and $\sigma=0.02$. On the two-qubit gates we
sample from a normal distribution with $\mu=0$ and $\sigma=0.005$. These are
chosen similarly to the errors in Ref. Niu  _et al._ (2022).

For the generator we use a layered circuit Ansatz such that
$G(\bm{\theta}_{g})=G_{L}(\bm{\theta}_{L})\cdots
G_{2}(\bm{\theta}_{2})G_{1}(\bm{\theta}_{1})$, where each layer,
$G_{l}(\bm{\theta}_{l})$, consists of an Euler rotation on each qubit followed
by nearest-neighbor couplings of all qubits using two-qubit gates. This means
that we have $3N$ single-qubit rotations in each layer, where $N$ is the
number of qubits. An example of a layer can be seen in Fig. 2, where the layer
is shown with cnot gates.

We apply the EQ-GAN to several different real quantum states. First, we
consider learning a simple two-qubit circuit with the same structure as the
generator, then we consider states generated by a VQE algorithm, and finally,
we consider completely random states.

###  II.1 Random two-qubit circuit

Figure 3: Circuit layout of a two-qubit EQ-GAN with a generator using the same
circuit layout as the real quantum state. Note that the generator and real
circuit constitute one layer of the circuit Ansatz shown in Fig. 2.

First, we consider a simple circuit with a one-layer layout identical to the
generator, but with randomly chosen predefined rotation angles,
$\bm{\theta}^{*}$. By choosing the layout of the generator to be identical to
the real data, we should, in principle, be able to generate the exact same
state. The full layout of the EQ-GAN can be seen in Fig. 3.

We train for 80 episodes and a batch size of 4, with learning rates of 0.01
for both networks. We sample ${10}^{4}\text{\,}$ different real parameters,
$\bm{\theta}^{*}$, and perform an EQ-GAN simulation on each sample for both a
perfect SWAP test and an adversarial SWAP test. The distribution of
fidelities, i.e.,
$|\bra{0}G(\bm{\theta}_{g})^{\dagger}R(\bm{\theta}^{*})\ket{0}|^{2}$, where
$G(\bm{\theta})\ket{0}$ is the generated state and $R(\bm{\theta}^{*})\ket{0}$
is the real state, can be seen in Fig. 4 where we also have plotted the
average fidelity of all the simulations. We observe that the perfect swap
performs better than the adversarial EQ-GAN on average. The perfect swap
simulations have a larger peak and density above 0.95 and a smaller tail.
However, the mode, i.e., the most frequent fidelity, of the adversarial
training is larger than that of the perfect swap. The fact that the mode of
the data is significantly larger than the fidelity is means that the data are
more spread out as we also observe. This is a result of the difficulty in
training adversarial networks. The perfect SWAP test requires no training
besides the generator and thus it is easier to train on average. On the other
hand, the adversarial EQ-GAN is more difficult to train as it requires both
training of the generator and training of the discriminator, yielding a lower
average fidelity but a large mode. This larger mode shows that the EQ-GAN has
the potential to be better trained, i.e., yielding a higher maximum fidelity,
compared with the perfect SWAP training. This difficulty in training GANs is
well documented for the classical versions Mescheder  _et al._ (2018).

Figure 4: Distribution of ${10}^{4}\text{\,}$ EQ-GAN simulations using the
simple circuit Ansatz shown in Fig. 3. The dashed lines indicate the average
fidelity of the simulations for each case.

Table 1 shows the rotation angles of the best EQ-GAN result. Despite a
relatively high fidelity (above 0.9999), we observe that the rotation angles
are far from identical to the real rotation angles, $\bm{\theta}^{*}$. Since
it is possible to rotate around the Bloch sphere using different approaches
and still end up at the same place, this is not an immediate problem. However,
this is a consequence of the fact that the fidelity measure, and thus the SWAP
test, only considers the amplitude of a given state and not the phases. Thus
we have not generated two identical states; rather we have just generated two
states with identical amplitudes. This should be kept in mind when working
with EQ-GANs, as it is typical for all simulations.

Table 1: Parameters used in the EQ-GAN simulation and overlap fidelity of the states returned by the circuits. $\bm{\theta}^{*}$ refer to the predefined rotation angles of the real circuit in Fig. 3. $\bm{\theta}_{p}$ refer to rotation angles generated using a perfect SWAP test, and $\bm{\theta}_{a}$ refer to rotation angles generated using an adversarial SWAP test as the discriminator. All angles are in radians. | Fidelity | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | $\theta_{4}$ | $\theta_{5}$ | $\theta_{6}$
---|---|---|---|---|---|---|---
$\bm{\theta}^{*}$ | 1 | 0.035 | 2.861 | 0.606 | 0.361 | 6.174 | 4.513
$\bm{\theta}_{p}$ | 0.95066 | 0.316 | 0.546 | 0.295 | 5.644 | 0.662 | 5.663
$\bm{\theta}_{a}$ | 0.99995 | 0.219 | 6.273 | 0.347 | 5.589 | 0.056 | 5.567

###  II.2 VQE learning

We now increase the complexity of our EQ-GAN model, both in the number of
qubits and in the number of layers of the generator, and by using real data.
As in the previous section, we keep the circuit architecture identical for
both the real data and the generator, meaning that in principle we should be
able to generate a circuit with unity fidelity. However, instead of learning
random states as in Ref. Benedetti  _et al._ (2019), we want to learn states
produced by a variational quantum eigensolver Peruzzo  _et al._ (2014);
McClean  _et al._ (2016); O’Malley  _et al._ (2016); Kandala  _et al._ (2017);
Cao  _et al._ (2019); Barkoutsos  _et al._ (2018); McCaskey  _et al._ (2019);
Gard  _et al._ (2020). In particular, we used results from an subspace-search
variational quantum eigensolver (SSVQE) Nakanishi  _et al._ (2019) simulation
used for approximating small benchmark molecules The reason we use this
particular variational quantum algorithm is that it can encode several states
in a single circuit, thus encoding more information in the circuit. The SSVQE
model has the same layered circuit Ansatz as the EQ-GAN, i.e., the one in Fig.
2, with the same number of layers. Simulating the SSVQE algorithm for the two
lowest eigenstates of a molecular Hamiltonian produces a set of parameters
$\bm{\theta}_{\text{VQE}}$. The real data are thus generated by the circuit
represented by $R(\bm{\theta}_{\text{VQE}})$. In other words, we are trying to
learn approximate eigenstates of a given Hamiltonian, without knowing the
actual Hamiltonian.

When the circuit we are trying to learn becomes bigger we need more data to
make the training more efficient. Therefore we generate more training data
from each VQE-eigenstate by sampling 100 sets of parameters from a normal
distribution with $\bm{\theta}_{\text{VQE}}$ as the mean value of a standard
deviation of 0.01. Other than that we keep the hyperparameters identical to
the ones used in Section II.1. We use the circuit Ansatz shown in Fig. 2.

In Fig. 5 we present results for simulations of H2. In Appendix B, we show
results for LiH and BeH2. The active space of H2 and LiH is encoded using four
qubits, while the active space of BeH2 is encoded using six qubits. We
simulate the molecules for different atomic bond lengths performing a
simulation for each $0.1\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. We plot the
energy of the VQE eigenstates and the energy of the generated states. We also
plot the infidelity,
$1-|\bra{0}G(\bm{\theta}_{g})^{\dagger}R(\bm{\theta}_{\text{VQE}})\ket{0}|^{2}$,
as a function of the bond length. Starting from the infidelity in Fig. 5(b),
we see that the adversarial training outperforms the perfect SWAP test, and in
general, we find an infidelity on the order of $10^{-2}$ equivalent to a
fidelity around 0.99. Despite this quite good infidelity, we find that the
energy of the generated states in Fig. 5(a) is nowhere near the SSVQE
energies. They are consistently above the SSVQE energies for both the ground
state and the first excited state, no matter which discriminator we use. This
is due to the lack of phase estimation of the EQ-GAN. As pointed out, the lack
of phase estimation is a limitation of the EQ-GAN. However, it is not likely
to be something that can be mended as that may imply a violation of the no-
cloning theorem. The phases of an eigenstate are important, and as the EQ-GAN
does not approximate these the resulting energies are poorly reproduced,
despite the overlap with the VQE-eigenstates. For LiH and BeH2 we find similar
results, see Appendix B.

Figure 5: Results of the EQ-GAN learning on VQE data for H2 with three layers
of the circuit Ansatz shown in Fig. 2. (a) Predicted energy of the EQ-GAN
algorithm (lines) and SSVQE (circles). Solid lines indicate that adversarial
training was used, while the dashed line indicate that the perfect SWAP test
was used. (b) Infidelity of the state generated by the EQ-GAN,
$\ket{\textup{GAN}}=G(\bm{\theta})\ket{0}$, with the state generated by the
VQE, $\ket{\textup{VQE}}=R(\bm{\theta}_{\textup{VQE}})\ket{0}$.

###  II.3 Random state learning

Finally, we consider how efficiently an EQ-GAN can learn random states. Thus
$R$ is now a completely random unitary, which is not known to the generator
circuit. This means that there is no guarantee that the generator will be able
to reproduce $R$ efficiently. Since we use the circuit Ansatz as in Fig. 2,
which only employs nearest-neighbor interactions, it is expected that the
fidelity will decrease as the number of qubits increases. Instead of just
using cnot gates as the entangling operation we also simulate the EQ-GAN using
imaginary swap ($i$swap) and cz gates and their parametrized equivalents
Rasmussen and Zinner (2022). Though it may seem like an advantage that the
parametrized two-qubit gates have more parameters, it is not. This is because
we have already saturated the number of parameters as discussed in Ref.
Rasmussen  _et al._ (2020).

We keep the training parameters as in the previous sections and sample 200
random states. For each random state, we sample 100 training states using the
same approach as in Section II.2, i.e., from a normal distribution with the
coefficients of the original state as the mean value and a standard deviation
of 0.01.

Figure 6: Average infidelity of the generated states for EQ-GAN when trying to
learn random states. The sizes of the generated states are as follows: (a) two
qubits, (b) four qubits, and (c) six qubits. Solid symbols indicate the fixed
version of the two-qubit gate, while open symbols indicate parametrized
versions of the two-qubit gate.

We present the average infidelity of the generated states in Fig. 6. For the
two-qubit case in Fig. 6(a), we observe an average infidelity of 0.1 for the
worst cases and below 0.01 for the best cases. The cz gate and its
parameterize version do not perform very well, and they actually perform worse
for increasing layers. This seems to be the case also for four and six qubits.
Thus we conclude that the cz gate is the least suited two-qubit gate for
generative learning. This is likely because the cz gate is diagonal in the
computational basis.

The $i$swap gate performs better than the cz gate, however, in the two-qubit
case, its performance is approximately constant at 0.04 as the number of
layers increases. The parametrized version performs better close to 0.01 on
average. The cnot gate also yields an infidelity around 0.01 on average for
more than two layers, indicating that the maximum average fidelity is likely
reached for this circuit Ansatz. For one layer we see a slight advantage of
the parametrized version of the cnot gate compared with the usually fixed
gate.

In the four- and six-qubit cases [Fig. 6(b) and (c)], we see an increased
minimum infidelity, since we now only have nearest-neighbor interactions. For
the four-qubit case, we find the lowest infidelity at 0.02, whereas the six-
qubit case does not perform better than 0.1. This increase in infidelity could
be mended by changing the circuit architecture; see e.g., Ref. Sim  _et al._
(2019) for a discussion of different circuit architectures. Again the cnot
gate performs best, but the infidelity is decreasing as the number of layers
is increasing. The fixed and parametrized cnot gates perform approximately
identically. For the $i$swap gate the parametrized version outperforms the
fixed version for one or two layers while they perform identically for more
layers. This indicates the advantages of using parametrized two-qubit gates
instead of just increasing the number of layers. This could be advantageous in
near-term devices Rasmussen and Zinner (2022).

##  III Conclusion

In this paper, we have considered the EQ-GAN and its ability to learn
different unitaries. We have shown that the EQ-GAN can perform better than the
SWAP test, but it is more difficult to train. The perfect SWAP test yields a
higher average fidelity, however, the EQ-GAN yields a higher maximum fidelity.
We also discuss that despite this high fidelity the trained parameters are not
necessarily the same as the real parameters used to define the state. This is
because the EQ-GAN lacks phase estimation of the generated states. We show
that this may be a problem, e.g., when trying to learn eigenstates of a
Hamiltonian. In such a case, even with high fidelity, the eigenvalues of the
generated states can be far from the target eigenvalues. Nonetheless, EQ-GANs
can still be useful if phase estimation is not required, for instance, if we
are looking for the probability of finding the amplitude of different basis
states in a given quantum states.

Finally, we discussed random state learning with a layered nearest-neighbor
circuit Ansatz. We find that the cnot gate performs best and that the
parametrized versions of two-qubit gates outperform the fixed versions for a
few layers. We also find that the fidelity decreases as the number of qubits
in the state we want to generate increases. This could be mended by changing
the circuit Ansatz to increase the qubit interactions beyond nearest neighbor.
The results could also be improved by tuning the hyperparameters of the model.
In general, hyperparameter tuning is a large issue in classical generative
models Salimans  _et al._ (2016), and it is therefore expected to be of
similar difficulty for quantum generative models.

Possible application of EQ-GANs could be to test how well various circuit
Ansätze can approximate sets of random states. Our results also show that EQ-
GANs could be useful for quantum state loading or distribution loading Zoufal
_et al._ (2019); Stamatopoulos  _et al._ (2020).

###### Acknowledgements.

The authors thank T. Bækkegaard for discussions on different aspects of the
work. This work is supported by the Danish Council for Independent Research.

##  Appendix A Cost function of EQ-GAN

A classical GAN consists of two neural networks, a generative network,
$G(\bm{\theta}_{g},\bm{z})$ and a discriminative network,
$D(\bm{\theta}_{d},\bm{z})$. The generator maps a vector, $\bm{z}$, sampled
randomly from a given distribution to a data example,
$G(\bm{\theta}_{g},\bm{z})$. We call this the fake data. The discriminator
takes a sample $\bm{x}$ and returns the probability,
$D(\bm{\theta}_{d},\bm{x})$, for the data being real or fake, i.e., generated
by the generator network.

This creates a minimax optimization problem, where we are alternating between
improving the discriminator’s ability to distinguish between real and fake
data and improving the generator’s ability to trick the discriminator. In
general, a classical GAN is trained by solving
$\min_{\bm{\theta}_{g}}\max_{\bm{\theta}_{d}}V(\bm{\theta}_{d}\bm{\theta}_{g})$,
where the cost function is given as

$\displaystyle V(\bm{\theta}_{d},\bm{\theta}_{g})=$
$\displaystyle\mathbb{E}_{x\sim p_{\text{data}}(x)}[\log
D(\bm{\theta}_{d},x)]$ (1) $\displaystyle+$ $\displaystyle\mathbb{E}_{z\sim
p_{0}(z)}\\{\log[1-D(\bm{\theta}_{d},G(\bm{\theta}_{g},z))]\\},$

where $\mathbb{E}_{x\sim p_{\text{data}}(x)}$ and $\mathbb{E}_{z\sim
p_{0}(z)}$ represent the expectation over the distributions $p_{\text{data}}$
and $p_{0}$ respectively. If the two neural networks approach the space of
arbitrary functions, then the global optimum exists and is uniquely determined
as $p_{g}(z)=p_{\text{data}}(z)$.

The first proposal of a QuGAN Lloyd and Weedbrook (2018); Dallaire-Demers and
Killoran (2018) defined the generative network as a quantum circuit,
$G(\bm{\theta}_{g})$, that outputs the state (here represented as a density
matrix) $\rho=G(\bm{\theta}_{g})\rho_{0}G^{\dagger}(\bm{\theta}_{g})$ from the
initial state $\rho_{0}$. The discriminator then takes a real state, $\sigma$,
or a fake state, $\rho$, as its input $\rho_{\text{in}}$, and performs a
positive operator-valued measurement (POVM) defined by $T$. The outcome of the
POVM determines the probability of the state being real,
$D(\bm{\theta}_{d},\rho_{\text{in}})=\tr[T\rho_{\text{in}}]$. This means that
the QuGAN is trained by solving the minimax problem

$\min_{\bm{\theta}_{g}}\max_{T}V(\bm{\theta}_{g},T)=\min_{\bm{\theta}_{g}}\max_{T}(\tr[T\sigma]-\tr[T\rho(\bm{\theta}_{g})]).$
(2)

However, this might not converge to the desired Nash equilibrium, leading to
oscillations of the generator and discriminator, as shown in Ref. Niu  _et
al._ (2022).

To fix this problem, the authors of Ref. Niu  _et al._ (2022) propose to let
the discriminator evaluate the true data, $\sigma$, and the fake data,
$\rho(\bm{\theta}_{g})$, simultaneously, as seen in Fig. 3. They define a
minimax cost function similar to the classical GAN by Eq. 1:

$V(\bm{\theta}_{g},\bm{\theta}_{d})=1-D(\bm{\theta}_{d},\rho(\bm{\theta}_{g})),$
(3)

where $D(\bm{\theta}_{d},\rho(\bm{\theta}_{g}))$ is the parametrized SWAP test
as seen in Fig. 1 and $\rho(\bm{\theta}_{g})$ is the generated state. Equation
3 defines the cost function used for the EQ-GAN considered in this paper.

##  Appendix B Additional results

Here we present some additional results. We present simulations for LiH with
one layer (Fig. 7) and BeH2 with three layers (Fig. 8).

Figure 7: Results of the EQ-GAN learning on VQE data for LiH with one layer.
(a) Predicted energy of the EQ-GAN algorithm (lines) and SSVQE (circles).
Solid lines indicate that the adversarial training was used, while the dashed
line indicate that the perfect SWAP test was used. (b) Infidelity of the state
generated by the EQ-GAN, $\ket{\textup{GAN}}=G(\bm{\theta})\ket{0}$, with the
state generated by the VQE,
$\ket{\textup{VQE}}=R(\bm{\theta}_{\textup{VQE}})\ket{0}$. Figure 8: Results
of the EQ-GAN learning on VQE data for BeH2 with three layers. (a) Predicted
energy of the EQ-GAN algorithm (lines) and SSVQE (circles). Solid lines
indicate that the adversarial training was used, while the dashed lines
indicate that the perfect SWAP test was used. (b) Infidelity of the state
generated by the EQ-GAN, $\ket{\textup{GAN}}=G(\bm{\theta})\ket{0}$, with the
state generated by the VQE,
$\ket{\textup{VQE}}=R(\bm{\theta}_{\textup{VQE}})\ket{0}$.

## References

  * Goodfellow  _et al._ (2014) I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, in  _Advances in Neural Information Processing Systems_, Vol. 27, edited by Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger (Curran Associates, Inc., 2014) pp. 2672–2680.
  * Karras  _et al._ (2019) T. Karras, S. Laine, and T. Aila, “A style-based generator architecture for generative adversarial networks,” (2019), arXiv:1812.04948 [cs.NE] .
  * Wang  _et al._ (2018) X. Wang, K. Yu, S. Wu, J. Gu, Y. Liu, C. Dong, C. C. Loy, Y. Qiao, and X. Tang, “ESRGAN: Enhanced super-resolution generative adversarial networks,” (2018), arXiv:1809.00219 [cs.CV] .
  * Mustafa  _et al._ (2019) M. Mustafa, D. Bard, W. Bhimji, Z. Lukic, R. Al-Rfou, and J. M. Kratochvil, Computational Astrophysics and Cosmology 6 (2019), 10.1186/s40668-019-0029-9.
  * Schawinski  _et al._ (2017) K. Schawinski, C. Zhang, H. Zhang, L. Fowler, and G. K. Santhanam, Monthly Notices of the Royal Astronomical Society: Letters 467 (2017), 10.1093/mnrasl/slx008.
  * Arjovsky  _et al._ (2017) M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein GAN,” (2017), arXiv:1701.07875 [stat.ML] .
  * Salimans  _et al._ (2016) T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, and X. Chen, “Improved techniques for training GANs,” (2016), arXiv:1606.03498 [cs.LG] .
  * Mescheder  _et al._ (2018) L. Mescheder, A. Geiger, and S. Nowozin, “Which training methods for GANs do actually converge?” (2018), arXiv:1801.04406 [cs.LG] .
  * Dallaire-Demers and Killoran (2018) P.-L. Dallaire-Demers and N. Killoran, Phys. Rev. A 98, 012324 (2018).
  * Lloyd and Weedbrook (2018) S. Lloyd and C. Weedbrook, Phys. Rev. Lett. 121, 040502 (2018).
  * Osborne and Rubinstein (1994) M. J. Osborne and A. Rubinstein,  _A Course in Game Theory_ (MIT Press, Cambridge, MA, 1994).
  * Niu  _et al._ (2022) M. Y. Niu, A. Zlokapa, M. Broughton, S. Boixo, M. Mohseni, V. Smelyanskyi, and H. Neven, Phys. Rev. Lett. 128, 220505 (2022).
  * Benedetti  _et al._ (2019) M. Benedetti, E. Grant, L. Wossnig, and S. Severini, New Journal of Physics 21, 043023 (2019).
  * Huang  _et al._ (2021) H.-L. Huang, Y. Du, M. Gong, Y. Zhao, Y. Wu, C. Wang, S. Li, F. Liang, J. Lin, Y. Xu, R. Yang, T. Liu, M.-H. Hsieh, H. Deng, H. Rong, C.-Z. Peng, C.-Y. Lu, Y.-A. Chen, D. Tao, X. Zhu, and J.-W. Pan, Phys. Rev. Applied 16, 024051 (2021).
  * Hu  _et al._ (2019) L. Hu, S.-H. Wu, W. Cai, Y. Ma, X. Mu, Y. Xu, H. Wang, Y. Song, D.-L. Deng, C.-L. Zou, and L. Sun, Science Advances 5, 10 (2019).
  * Zoufal  _et al._ (2019) C. Zoufal, A. Lucchi, and S. Woerner, npj Quantum Information 5, 103 (2019).
  * Situ  _et al._ (2020) H. Situ, Z. He, Y. Wang, L. Li, and S. Zheng, Information Sciences 538, 193–208 (2020).
  * Romero and Aspuru-Guzik (2019) J. Romero and A. Aspuru-Guzik, “Variational quantum generators: Generative adversarial quantum machine learning for continuous distributions,” (2019), arXiv:1901.00848 [quant-ph] .
  * Zeng  _et al._ (2019) J. Zeng, Y. Wu, J.-G. Liu, L. Wang, and J. Hu, Phys. Rev. A 99, 052306 (2019).
  * Chang  _et al._ (2021) S. Y. Chang, S. Herbert, S. Vallecorsa, E. F. Combarro, and R. Duncan, EPJ Web of Conferences 251, 03050 (2021).
  * Stein  _et al._ (2021) S. A. Stein, B. Baheri, D. Chen, Y. Mao, Q. Guan, A. Li, B. Fang, and S. Xu, “QuGAN: A generative adversarial network through quantum states,” (2021), arXiv:2010.09036 [quant-ph] .
  * Carleo and Troyer (2017) G. Carleo and M. Troyer, Science 355, 602 (2017), https://www.science.org/doi/pdf/10.1126/science.aag2302 .
  * Carleo  _et al._ (2017) G. Carleo, L. Cevolani, L. Sanchez-Palencia, and M. Holzmann, Phys. Rev. X 7, 031026 (2017).
  * Hibat-Allah  _et al._ (2020) M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Phys. Rev. Research 2, 023358 (2020).
  * Che  _et al._ (2022) Y. Che, C. Gneiting, and F. Nori, Phys. Rev. B 105, 214205 (2022).
  * Gupta and Zia (2002) S. Gupta and R. K. P. Zia, “Quantum neural networks,” (2002), arXiv:quant-ph/0201144 [quant-ph] .
  * Farhi and Neven (2018) E. Farhi and H. Neven, “Classification with quantum neural networks on near term processors,” (2018), arXiv:1802.06002 [quant-ph] .
  * Beer  _et al._ (2020) K. Beer, D. Bondarenko, T. Farrelly, T. J. Osborne, R. Salzmann, D. Scheiermann, and R. Wolf, Nature Communications 11 (2020), 10.1038/s41467-020-14454-2.
  * Buhrman  _et al._ (2001) H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Phys. Rev. Lett. 87, 167902 (2001).
  * Garcia-Escartin and Chamorro-Posada (2013) J. C. Garcia-Escartin and P. Chamorro-Posada, Phys. Rev. A 87, 052330 (2013).
  * Broughton  _et al._ (2021) M. Broughton, G. Verdon, T. McCourt, A. J. Martinez, J. H. Yoo, S. V. Isakov, P. Massey, R. Halavati, M. Y. Niu, A. Zlokapa, E. Peters, O. Lockwood, A. Skolik, S. Jerbi, V. Dunjko, M. Leib, M. Streif, D. V. Dollen, H. Chen, S. Cao, R. Wiersema, H.-Y. Huang, J. R. McClean, R. Babbush, S. Boixo, D. Bacon, A. K. Ho, H. Neven, and M. Mohseni, “TensorFlow Quantum: A software framework for quantum machine learning,” (2021), arXiv:2003.02989 [quant-ph] .
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferova, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Chemical Reviews 119, 10856 (2019).
  * Barkoutsos  _et al._ (2018) P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo, S. Filipp, and I. Tavernelli, Phys. Rev. A 98, 022322 (2018).
  * McCaskey  _et al._ (2019) A. J. McCaskey, Z. P. Parks, J. Jakowski, S. V. Moore, T. D. Morris, T. S. Humble, and R. C. Pooser, npj Quantum Information 5, 99 (2019).
  * Gard  _et al._ (2020) B. T. Gard, L. Zhu, G. S. Barron, N. J. Mayhall, S. E. Economou, and E. Barnes, npj Quantum Information 6, 10 (2020).
  * Nakanishi  _et al._ (2019) K. M. Nakanishi, K. Mitarai, and K. Fujii, Phys. Rev. Research 1, 033062 (2019).
  * Rasmussen and Zinner (2022) S. E. Rasmussen and N. T. Zinner, “Parameterized two-qubit gates for enhanced variational quantum eigensolver,” (2022), arXiv:2203.04978 [quant-ph] .
  * Rasmussen  _et al._ (2020) S. E. Rasmussen, N. J. S. Loft, T. Bækkegaard, M. Kues, and N. T. Zinner, Advanced Quantum Technologies 3, 2000063 (2020).
  * Sim  _et al._ (2019) S. Sim, P. D. Johnson, and A. Aspuru-Guzik, Advanced Quantum Technologies 2, 1900070 (2019).
  * Stamatopoulos  _et al._ (2020) N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, Quantum 4, 291 (2020).


# Discontinuous transition to shear flow turbulence

Y. Zhuang1 B. Yang1 V. Mukund1 E. Marensi2 & B. Hof1

###### Abstract

Depending on the type of flow the transition to turbulence can take one of two
forms, either turbulence arises from a sequence of instabilities, or from the
spatial proliferation of transiently chaotic domains, a process analogous to
directed percolation. Both scenarios are inherently continuous and hence the
transformation from ordered laminar to fully turbulent fluid motion is only
accomplished gradually with flow speed. Here we show that these established
transition types do not account for the more general setting of shear flows
subject to body forces. By attenuating spatial coupling and energy transfer,
spatio-temporal intermittency is suppressed and with forcing amplitude the
transition becomes increasingly sharp and eventually discontinuous. We argue
that the suppression of the continuous range and the approach towards a first
order, discontinuous scenario applies to a wide range of situations where in
addition to shear, flows are subject to e.g. gravitational, centrifugal or
electromagnetic forces.

Institute of Science and Technology Austria, 3400 Klosterneuburg, Austria

The University of Sheffield, Sheffield, UK

One of the earliest model for the onset of turbulence goes back to Landau[17]
who envisioned a sequence of instabilities that gives rise to a large number
of incommensurate modes and a stepwise increase of the flow’s complexity. This
view was later challenged and refined by Ruelle and Takens[26] who showed that
chaotic dynamics already arises after a small number of bifurcations. Such
universal routes[7] to chaos are typically observed in flows dominated by body
forces[9, 19, 31] (e.g. buoyancy or centrifugal forces) and subject to a
linear instability of the base state.

Conversely, in shear flows including those in pipes, channels and boundary
layers, the transition is fundamentally different. Turbulence tends to result
from perturbations of finite amplitude and the resulting dynamics are
immediately high dimensional. Specifically the turbulent state does not
develop from a linear instability, but from nonlinear, unstable solutions that
are unrelated and dynamically disconnected from the laminar base flow[14, 6,
13, 10, 23]. Despite the abrupt switch between laminar and turbulent dynamics,
the transition is nevertheless continuous. Starting from fully turbulent flow,
turbulence does not disappear abruptly with decreasing Reynolds number
($Re=UD/\nu$, where for pipe flow $U$ is chosen as the bulk velocity, $D$ the
pipe diameter and $\nu$ the fluid’s kinematic viscosity). Instead turbulence
ceases to be space filling and the fraction of the flow that is turbulent
decreases gradually until it eventually reaches zero at the transition
threshold .

The co-existence of laminar and turbulent regions and the associated spatio-
temporally intermittent dynamics arise from spatial coupling and energy
transfer between laminar and turbulent regions. In pipe flow, turbulent
structures of limited size, puffs, prevail by extracting energy from the
upstream laminar flow[32, 33]. The high speed laminar fluid close to the
centreline impinges onto the adjacent puff and gives rise to a peak in
turbulent production at the puff’s upstream interface[12]. Conversely, if this
energy influx is suppressed, by distorting the shape of the upstream laminar
profile from a parabolic to a plug shape (using a ’plug forcing’ scheme[12]),
the production peak disappears and the puff decays.

Even without manipulation of the upstream flow, puffs eventually decay. Yet,
prior to decay they can split and give rise to new puffs[1]. This competition
between decay and proliferation determines the aforementioned transition
threshold and sets the stage for a directed percolation (i.e. continuous)
phase transition. This however is not the end of the transition process. Puffs
cannot merge to form larger continuous structures and the flow remains spatio
temporally intermittent. The delocalization of puffs and the eventual rise of
fully turbulent flow[4] can be understood by viewing pipe flow as an
excitable-bistable medium [2, 3] or in a related approach as actuator-
inhibitor dynamics[28, 8, 35]. Just like the preceding DP phase transition,
the subsequent emergence of fully turbulent flow[4] is also a continuous
process, i.e. it occurs gradually with Reynolds number.

The two transition types described above are respectively referred to as
super- and subcritical transition. It is important to emphasize that this
nomenclature takes the perspective of the laminar state, not that of
turbulence. In the former case this perspective is indeed the appropriate one
because the transition indeed follows from a supercritical bifurcation of the
laminar state. In the latter case however, the stability of the laminar flow
and a possible bifurcation at higher Re, which the term subcritical alludes
to, are entirely irrelevant for the formation of turbulence[23]. As outlined
above, the relevant transition threshold is set by the proliferation of
turbulence[1] and instead of the start of a bifurcation sequence it demarcates
a phase transition[27, 28, 18, 5, 11].

The supercritical transition scenario covers situations where flows are
destabilized by body forces, whereas the subcritical one prevails in shear
flows. However, most practical situations are more complicated and in addition
to shear, flows are equally subject to body forces. Examples range from basic
engineering applications such as heating pipes and flows around bends and
corners, all the way to geo- and astrophysical flows including the earth’s
boundary layer or flows of plasmas and ionized gases subject to magnetic
fields. Naively one would expect that for such more complicated scenarios the
transition is either dominated by shear or by a sequence of instabilities that
may arise from the additional force(s) and that in either case the result is a
continuous transition. As we will show in the following, body forces tend to
diminish spatial interactions and in doing so they severely suppress laminar
turbulence co-existence. Consequently the spatio-temporally intermittent
regime, believed intrinsic to the onset of shear flow turbulence, can be fully
eradicated giving way to a much simpler discontinuous scenario and a
transition directly from fully turbulent to fully laminar.

We start our investigation with flows through curved pipes and in this case
curvature is well known to delay transition[36, 30]. More precisely we perform
experiments in a helical pipe which is composed of a semi-flexible tube that
is coiled around a cylinder of radius A=22.5cm, while the pipe has an inner
radius of R=4mm, resulting in a radius ratio R/A=0.0173 and a pitch of the
resulting helix of 12 mm or 3R per coil. The tube has a length of L=98 m and a
dimensionless length L/R= 24 000. The addition of flow visualization particles
allows the discrimination between laminar and turbulent regions and velocity
fields were measured using particle image velocimetry (PIV). While pipe
curvature generally suppresses turbulence, for radius ratios $\gtrsim$ 0.25 a
linear instability[30] of the laminar flow occurs prior to the onset of
turbulence[15, 21, 25]. For our chosen radius ratio of 0.0173 this linear
instability does not play any role. Hence the transition remains ’subcritical’
and occurs in response to finite amplitude perturbations. In the present case
the flow was perturbed directly at the inlet: prior to entering the curved
section the fluid passes through a 300 R long straight tube of identical
diameter. At the beginning of the straight pipe a small pin protrudes into the
tube and assures fully turbulent inflow conditions provided that $Re\gtrsim
2800$.

In agreement with these earlier observations we find that curvature delays the
transitions. The turbulent inflow relaminarizes rapidly, within the first few
turns of the helical pipe, provided that $Re\lesssim 3400$. For increasing Re
turbulence persists for longer times and hence survives to further downstream
distances (note that turbulence advects downstream at a speed close to the
bulk velocity). A recent computational study[25] of transition in curved pipes
observed that in this geometry the turbulent kinetic energy peak at the
interface, and hence the signature of the spatial coupling between laminar and
turbulent domains, is strongly suppressed and that puff proliferation
(splitting) is absent. In addition to the absence of strong fronts (red signal
in Fig. 1c) not only do we not find any puff splittings but puffs were
entirely absent for the parameters used in experiments. Although short
localized patches can be triggered they are not metastable and either decay or
expand (depending on Re). Fig. 1a shows the fraction of the flow that is still
turbulent at five different downstream locations as a function of the Reynolds
number. These turbulent fraction curves are S-shaped and with downstream
location, they shift to higher Re. As illustrated by the flow visualization
images (Fig. 1b), turbulence continues to decay with downstream distance.
Notably the turbulent fraction curves become steeper, indicating that the
spatio-temporal intermittent regime is suppressed and that the transition
becomes increasingly sharp. Despite the length of more than 24 000 R the
turbulent fraction has not reached a statistical steady state at the end of
the pipe and apart from the trend towards discontinuity it is hence not
possible in the given set up to determine if the transition actually is
discontinuous in this case.

We next switch back to straight pipes and instead of curvature and the
resulting centrifugal force, we consider the effect of buoyancy on transition.
More specifically we carry out direct numerical simulations of a heated
vertical pipe (see reference[20] The flow is driven by a pressure gradient in
the upward direction and at the same time the fluid adjacent to the heated
wall is of lower density and experiences an additional upward driving. Like
curvature also buoyancy delays the transition to turbulence, a circumstance
that can significantly reduce the heat transfer in cooling pipes. We
initialize the flow in a fully turbulent state and (like in the curved pipe
case) monitor the turbulent fraction as a function of time. In Fig.1 c we show
the turbulent fraction encountered at three different time instances as a
function of Re. Also in this case the turbulent fraction curves are initially
S-shaped and steepen as time proceeds, qualitatively resembling the
phenomenology in curved pipes. In this case however, owing to the periodic
boundary conditions in the simulations, the observation times are not limited
by the physical length of the pipe. At each Re simulations were continued
until eventually the data approached a statistical steady state (blue data
points in Fig. 1c), often requiring more than $10^{4}$ advective time units
($R/(U_{cl})$). The data attests that albeit spatio temporal intermittency may
persists for long times, turbulent fractions $\lesssim 0.8$ are not sustained
and the transition is discontinuous. It is noteworthy that also for the heated
pipe puffs are absent.

The two body forces considered above affect the flow in very different ways,
while curvature gives rise to a centrifugal force and distorts the profile
radially, buoyancy accelerates the flow in the axial direction. Surprisingly
the impact on the transition to turbulence is qualitatively the same.

To better understand why in both cases spatio temporal intermittency is
suppressed and the transition becomes increasingly abrupt, we compare the
respective laminar and turbulent velocity profile and recall that STI arises
from spatial coupling and requires energy exchange between laminar and
turbulent regions. As depicted in Fig. 2a in the absence of body forces the
time averaged turbulent velocity profile is plug like and has a far smaller
peak velocity compared to the parabolic laminar flow. This large velocity
difference facilitates spatial coupling and energy flux from laminar to
turbulent regions. The situation is markedly different in the curved pipe case
(Fig. 2b). The centrifugal force affects both the laminar and the turbulent
flow and overall the effect of the body force on the profile exceeds the
profile deformation caused by the turbulent eddies. Hence the body force
effectively reduces the difference between the laminar and the turbulent
profile and in doing so it limits the spatial energy transfer between the two
states. As a consequence spatio temporal intermittency is suppressed. The
transition is delayed until at larger Re the locally available shear suffices
to sustain turbulence throughout the pipe not depending on spatial energy
flux. Equally in case of the heated pipe buoyancy affects both, the laminar
and turbulent profile shape. Again the profile distortion caused by buoyancy
far exceeds the distortion due to turbulence and hence also in this case the
laminar and the turbulent profiles show much smaller differences than for
ordinary pipe flow. The same argument applies, lacking spatial coupling and
energy transfer, spatio temporal intermittency is suppressed, rendering the
transition discontinuous.

This line of argument should equally apply to other body forces and hence the
suppression of the continuous range and the approach to a discontinuous
transition may be encountered in a broad range of situations. To probe this we
carried out simulations of two generic body forces that affect the streamwise
velocity profile in different ways. The first has the tendency to make the
profile more plug like while the second tends to preserve a parabolic profile
shape. The ’plug forcing’ has been used in a number of previous studies’ and
is known to suppress turbulence[12, 16]. The flattening of the velocity
profile is a property shared with magnetohydrodynamic pipe flow, where a
transverse magnetic field induces a Lorentz force that decelerates the central
flow and accelerates the near wall fluid. Comparable average velocity profiles
are equally found in pipe flows of shear thinning fluids. The ’parabolic
forcing’ scheme has the opposite effect on turbulent flow and accelerates the
flow at the pipe centre and tends to decelerate the near wall fluid. More
precisely the amplitude of the force is proportional to the local deviation
from the parabolic profile, with the proportionality constant set to
$\alpha=0.2$. Hence, unlike the other forcing schemes, this body force
actively reduces the difference between the laminar and the turbulent profile
and it acts locally and instantaneously. As shown in Fig. 2d,e also the plug
and parabolic forcing severely reduce the difference between the forced
laminar and the forced turbulent profile (compared to Fig. 2a). As a result
the energy peak at the upstream laminar turbulent interfaces is significantly
reduced and strong fronts are absent as shown in Fig. 2f.

In order to quantify the energy input from laminar flow we compute the energy
flux between laminar and turbulent regions for all the body forces
investigated(Fig. 2g). A positive value denotes a net energy flux into the
turbulent region (see Supplemental information for definition). This quantity
is drastically reduced in the presence of the body forces tested. In the
heated pipe it is decreased by more than a factor of five when compared to the
unforced case and for the parabolic force the flux has almost vanished. For
the plug force the flux has become negative, which in this case is caused by a
structural change of turbulence, which is confined to the strong near wall
shear layers, a situation that is equally observed in MHD pipe flows. Either
way the kinetic energy provided by the adjacent laminar flow is strongly
reduced in all cases. We would therefore expect that as in the case of the
heated pipe also the plug and parabolic forcing suppress STI and cause a
sharper if not discontinuous transition.

In order to be able to compare the abruptness of the transitions encountered
for the various body forces, we display in Fig. 3a the turbulent fraction as a
function of the reduced Reynolds number, $\epsilon=(Re-Re_{c})/Re_{c}$, where
$Re_{c}$ denotes the respective critical point. Flows were initially fully
turbulent and the experiment and simulations were then continued until a
statistical steady state was approached. In addition measurements were carried
out where the initial flow field was composed of laminar and turbulent regions
(typically TF=0.5) to detect for hysteresis and metastability. Compared to
unforced pipe flow all four body forces display a sharp transition and the
turbulent fraction increases from zero to one across narrow parameter range as
shown in Fig. 3a. In fact in all cases the change from $TF=0$ to $TF=1$ is
accomplished within less than $\epsilon=(Re-Re_{c})/Re_{c}=0.03$ above
critical, compared to $\epsilon\approx 0.5$ for ordinary pipe flow . For the
curved and heated pipe cases values $0<TF<1$ could be measured within this
narrow Re regime. For the plug forcing no intermittency could be detected and
to within $\epsilon\approx 0.002$ the transition was direct from laminar to
fully turbulent.

We will next focus on the parabolic forcing. In this case the increase in TF
from zero to one, occurs within a Reynolds number interval of less than
$\Delta Re=0.25$ corresponding to $\epsilon=0.0001$ or $0.01\%$ of the
critical value. The sharp switch between states and the discontinuous nature
of this transition is further underlined by strong hysteresis as illustrated
in Fig.3b. Coming from high Reynolds numbers fully turbulent flow persists far
below the critical point (open triangles in Fig.3b). In these cases the fully
turbulent flow is found to be metastable. Once a laminar gap opens up (see
Figure 3c for an example at Re=2650), the laminar nucleus grows and just like
a seed crystal in a super cooled liquid, it completely eliminates the
metastable phase. The transition hence shows all the characteristics expected
at a first order phase transition.

From its historical beginnings[24, 22, 29] the investigation of the transition
to turbulence has revolved around instabilities, and bifurcation sequences
were at the core of the conceptual approaches and theories[17, 26]. By
exploring the effect of body forces on the onset of turbulence in shear flows
we underline that the central question is not one of instability, but of
sustainability[34] of turbulence and of the associated phase transition[1].
While this framework has been applied in a variety of recent studies[28, 18,
5, 11], in all cases so far the phase transition was not only continuous but
also it merely manifested the initial step in a complex spatio-temporal
process. The transformation to fully turbulent flow appears as a separate
problem[4, 3] and is only established after an additional transition. In this
sequence the initial phase transition may appear from practical considerations
the less relevant step, as the main drag increase occurs during the later
transformation only. Our study demonstrates that body forces can greatly
simplify the transition scenario and eradicate the spatio-temporal complexity.
What remains is a discontinuous phase transition which, given the finite
amplitude separation between the turbulent and the laminar state is the
simplest and most natural scenario and in that respect can be regarded as the
baseline of the transition to turbulence.

From a practical point of view, body forces acting on shear flows are the rule
rather than the exception. Our study demonstrates that in these situations the
continuous regime, which based on the standard subcritical transition scenario
is expected to dominate flows, will often be significantly reduced and may
approach the discontinuous case.

## References

  * [1] Kerstin Avila, David Moxey, Alberto de Lozar, Marc Avila, Dwight Barkley, and Björn Hof.  The onset of turbulence in pipe flow.  Science, 333(6039):192–196, 2011.
  * [2] Dwight Barkley.  Simplifying the complexity of pipe flow.  Phys. Rev. E, 84(1):016309, 2011.
  * [3] Dwight Barkley.  Theoretical perspective on the route to turbulence in a pipe.  J. Fluid Mech., 803:P1, 2016.
  * [4] Dwight Barkley, Baofang Song, Vasudevan Mukund, Grégoire Lemoult, Marc Avila, and Björn Hof.  The rise of fully turbulent flow.  Nature, 526(7574):550–553, 2015.
  * [5] Matthew Chantry, Laurette S Tuckerman, and Dwight Barkley.  Universal continuous transition to turbulence in a planar shear flow.  Journal of Fluid Mechanics, 824, 2017.
  * [6] Bruno Eckhardt, Tobias M Schneider, Bjorn Hof, and Jerry Westerweel.  Turbulence transition in pipe flow.  Annu. Rev. Fluid Mech., 39:447–468, 2007.
  * [7] Mitchell J Feigenbaum.  Quantitative universality for a class of nonlinear transformations.  Journal of statistical physics, 19(1):25–52, 1978.
  * [8] Nigel Goldenfeld and Hong-Yan Shih.  Turbulence as a problem in non-equilibrium statistical mechanics.  Journal of Statistical Physics, 167(3-4):575–594, 2017.
  * [9] J. P. Gollub and Harry L. Swinney.  Onset of turbulence in a rotating fluid.  Phys. Rev. Lett., 35:927–930, Oct 1975.
  * [10] Michael D. Graham and Daniel Floryan.  Exact coherent states and the nonlinear dynamics of wall-bounded turbulent flows.  Annu. Rev. Fluid Mech., 53(1):227–253, 2021.
  * [11] B. Hof.  Directed percolation and the transition to turbulence.  Nat. Rev. Phys., 5:62–72, 2023.
  * [12] Björn Hof, Alberto De Lozar, Marc Avila, Xiaoyun Tu, and Tobias M Schneider.  Eliminating turbulence in spatially intermittent flows.  Science, 327(5972):1491–1494, 2010.
  * [13] Genta Kawahara, Markus Uhlmann, and Lennaert Van Veen.  The significance of simple invariant solutions in turbulent flows.  Annu. Rev. Fluid Mech., 44:203–225, 2012.
  * [14] RR Kerswell.  Recent progress in understanding the transition to turbulence in a pipe.  Nonlinearity, 18(6):R17, 2005.
  * [15] Jakob Kühnen, P Braunshier, M Schwegel, HC Kuhlmann, and Björn Hof.  Subcritical versus supercritical transition to turbulence in curved pipes.  Journal of Fluid Mechanics, 770, 2015.
  * [16] Jakob Kühnen, Baofang Song, Davide Scarselli, Nazmi Burak Budanur, Michael Riedl, Ashley P Willis, Marc Avila, and Björn Hof.  Destabilizing turbulence in pipe flow.  Nature Physics, 14(4):386–390, 2018.
  * [17] Lev D Landau.  On the problem of turbulence.  In Dokl. Akad. Nauk USSR, volume 44, page 311, 1944.
  * [18] Grégoire Lemoult, Liang Shi, Kerstin Avila, Shreyas V Jalikop, Marc Avila, and Björn Hof.  Directed percolation phase transition to sustained turbulence in Couette flow.  Nature Physics, 12(3):254, 2016.
  * [19] Albert Libchaber and Jean Maurer.  Une expérience de rayleigh-bénard de géométrie réduite; multiplication, accrochage et démultiplication de fréquences.  Le Journal de Physique Colloques, 41(C3):C3–51, 1980.
  * [20] Elena Marensi, Shuisheng He, and Ashley P Willis.  Suppression of turbulence and travelling waves in a vertical heated pipe.  Journal of Fluid Mechanics, 919:A17, 2021.
  * [21] Azad Noorani and Philipp Schlatter.  Evidence of sublaminar drag naturally occurring in a curved pipe.  Physics of Fluids, 27(3), 2015.
  * [22] William M’F Orr.  The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. part ii: A viscous liquid.  In Proceedings of the Royal Irish Academy. Section A: Mathematical and Physical Sciences, volume 27, pages 69–138. JSTOR, 1907.
  * [23] Chaitanya S Paranjape, Gökhan Yalnız, Yohann Duguet, Nazmi Burak Budanur, and Björn Hof.  Direct path from turbulence to time-periodic solutions.  arXiv preprint arXiv:2306.05098, 2023.
  * [24] Rayleigh.  On the stability or instability of certain fluid motions, ii.  Proceedings of the London Mathematical Society, 1(1):67–75, 1887.
  * [25] Enrico Rinaldi, Jacopo Canton, and Philipp Schlatter.  The vanishing of strong turbulent fronts in bent pipes.  Journal of Fluid Mechanics, 866:487–502, 2019.
  * [26] David Ruelle and Floris Takens.  On the nature of turbulence.  Les rencontres physiciens-mathématiciens de Strasbourg-RCP25, 12:1–44, 1971.
  * [27] L. Shi, M. Avila, and B. Hof.  Scale Invariance at the Onset of Turbulence in Couette Flow.  Phys. Rev. Lett., 110(20):204502, 2013.
  * [28] H.-Y. Shih, T.-L. Hsieh, and N. Goldenfeld.  Ecological collapse and the emergence of travelling waves at the onset of shear turbulence.  Nat. Phys., 12(3):245–248, 2015.
  * [29] Arnold Sommerfeld.  Ein beitrag zur hydrodynamischen erklaerung der turbulenten fluessigkeitsbewegung.  Atti Congr. Int. Math. 4th Rome, 1908.
  * [30] KR Sreenivasan and PJ Strykowski.  Stabilization effects in flow through helically coiled pipes.  Experiments in Fluids, 1(1):31–36, 1983.
  * [31] Harry L Swinney and Jerry P Gollub.  The transition to turbulence.  Physics today, 31(8):41–49, 1978.
  * [32] Casimir van Doorne.  Stereoscopic PIV on transition in pipe flow.  PhD thesis, TU Delft, 2004.
  * [33] Casimir WH van Doorne and Jerry Westerweel.  The flow structure of a puff.  Phil. Trans. R. Soc. A, 367(1888):489–507, 2009.
  * [34] Fabian Waleffe.  Transition in shear flows. nonlinear normality versus non‐normal linearity.  Physics of Fluids, 7, 1995.
  * [35] X. Wang, H.-Y. Shih, and N. Goldenfeld.  Stochastic Model for Quasi-One-Dimensional Transitional Turbulence with Streamwise Shear Interactions.  Phys. Rev. Lett., 129(3):034501, 2022.
  * [36] CM White.  Streamline flow through curved pipes.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 123(792):645–663, 1929.

Figure 1: Approach towards a discontinuous transition. a Reynolds number
dependence of the turbulent fraction (TF) in a curved (helical) pipe in
experiments. Curves show the turbulent fractions measured at different
downstream locations. The turbulent fraction continues to adjust with
downstream location and has not reached the statistical steady state at the
end of the pipe (24 000 R). The adjustment from laminar (TF=0) to fully
turbulent (TF=1) flow occurs across a decreasing Re interval. b An initially
fully turbulent flow gradually laminarizes. Turbulent transients persist over
the first 8000 pipe radii until eventually the flow fully relaminarizes. c
Turbulent fraction as a function of Re, in a vertically heated pipe in direct
numerical simulations, Also in this case the turbulent fraction initially
adjusts and reaches a statistically steady state after approximately 10 000
advective time units. Turbulent fractions below 0.8 cannot be sustained and
the transition is hence discontinuous.  Figure 2: Comparison of laminar and
turbulent profiles. a (unforced) pipe flow, b curved pipe, c heated pipe, d
plug forcing, e parabolic forcing. f fluctuations of the wall normal kinetic
energy (q). The large energy peak for the unforced pipe case marks the
location of the strong upstream front. This strong upstream front is evident
for the unforced pipe case, but absent for all the body forces studied. g
Energy flux across the laminar turbulent interfaces. Figure 3: Discontinuity,
hysteresis and metastability. Panel a shows the turbulent fraction as a
function of the reduced Reynolds number $\epsilon$. Compared to ordinary
(unforced) pipe flow (black) the transition in the presence of the various
body forces is sharp. The inset highlights the width of the laminar turbulent
co-existence regimes (i.e. the $\epsilon$ interval in which $0<TF<1$). As
shown in panel b for the parabolic forcing the transition is hysteretic. Full
symbols show the turbulent fraction the flow settles down to, if initialized
with a turbulent fraction TF=0.5. Additionaly runs were carried out starting
from fully turbulent flow (open symbols) and in this case fully turbulent flow
persists below the critical point ($Re_{c}=2779$. As expected at a first order
phase transition, below the critical point turbulence is metastable and after
nucleation of a laminar gap (see panel c), the stable laminar state invades
turbulence and the flow fully laminarizes.


# Towards modified bimetric theories within non-product spectral geometry

Arkadiusz Bochniak  Institute of Theoretical Physics, Jagiellonian University,
prof. Stanisława Łojasiewicza 11, 30-348 Kraków, Poland.

###### Abstract

We discuss class of doubled geometry models with diagonal metrics. Based on
the analysis of known examples we formulate a hypothesis that supports
treating them as modified bimetric gravity theories. Certain steps towards the
generic case are then performed.

##  I Introduction

The description of gravity in terms of geometric objects is the cornerstone of
Einstein’s General Relativity and leads to an intriguing possibility of
geometrizing all of the fundamental interactions. One of the existing proposal
is based on the noncommutative geometry Co94  \- a framework that puts on
equal footing both the metric structure of manifolds, Yang-Mills-type theories
and also the Higgs mechanism. The spectral description of manifolds Co93  can
be generalized into other than classical geometries like discrete spaces and
their products with manifolds. The latter one leads to the definition of the
so-called almost-commutative geometries that were successfully applied to the
description of gauge theories Co96 ; Co96a ; ChaCoMa07 . Appropriate choice of
the finite space allows e.g. for the formulation of the noncommutative
Standard Model of Particle Physics. In this case the finite geometry is build
on the matrix algebra $\mathbb{C}\oplus\mathbb{H}\oplus M_{3}(\mathbb{C})$
whose choice is dictated by the gauge group of the model vSbook .

Yet another model of this type, but much more simpler, is the one studied by
Connes and Lott CoLo91 , where the finite algebra is just
$\mathbb{C}\oplus\mathbb{C}$ and corresponds to the two points. In this case,
the product space can be thought of as $M\times\mathbb{Z}_{2}$, that is, we
have two copies of the same manifold. One can further generalize this geometry
and can allow for two distinct metrics on these two sheets sitarz2019 . Such a
doubled geometry is beyond the usual almost-commutative framework and
therefore is of a non-product type. Since the spectral action principle
applied to a single copy produces the Hilbert-Einstein action, the natural
question of the form of an action functional for this non-product type of
geometries arises. The answer for generic choices of metrics is not known yet,
but in the case of the Friedmann-Lemaître-Robertson-Walker (FRLW) type of
Euclidean metrics this was done analitycally sitarz2019 ; BS2021 , and the
stability of certain solutions was also analysed BS2021 . It was demonstrated
therein that the interaction between the metrics resembles features
characteristic to bimetric gravity models HaRo12 ; AKMS13 . Despite numerous
similarities, certain significant differences are also present. In particular,
the interaction potential for bimetric model is a polynomial one, while for
the two-sheeted model it is a rational function. Further similarities and
differences for generic metrics were recently analysed in BoSi22 , where yet
another interpretation of this model in terms of interacting branes was
proposed.

In this note we discuss yet another class of models beyond the FLRW framework.
We illustrate the generally claimed features on the simplified example - the
so-called Hopf model. In this case the interaction potential has nontrivial
logarithmic terms but it still possesses bimetric gravity characteristics.
Finally, we make a general comment on the doubled geometry models which may
allow for its future numerical studies.

##  II The generic diagonal model

A framework of spectral geometry, allowing for an equivalent description of
geometric objects in terms of algebraic data, originates from the observation
that the geometry of a compact spin Riemannian manifold $M$ can be encoded in
the collection of data $(C^{\infty}(M),L^{2}(M),D_{M})$ Co93 , where
$L^{2}(M)$ is the Hilbert space of square-integrable spinors, and $D_{M}$ is
the Dirac operator, which can be written locally (with the use of the spin
connection $\omega$) as $i\gamma^{\mu}(\partial_{\mu}+\omega_{\mu})$.

This system of data is a prerequisite for the notion of a spectral triple, a
set $(A,H,D)$ consisting of a unital $\ast$-algebra $A$ represented in a
faithful way on a Hilbert space $H$, on which the (possibly unbounded) densely
defined self-adjoint operator $D$ acts. In the generic case it is assumed that
the commutators $[D,a]$, $a\in A$, are well-defined and can be (uniquely)
extended to an element from $B(H)$, bounded operators on $H$. Furthermore, the
resolvent of $D$ has to be compact. Several further comptability conditions
are imposed for certain applications Li18 ; vSbook .

In addition to the aforementioned canonical spectral triple associated to a
manifold $M$, the finite dimensional ones are well-understood Kra98 ; PaSi98 .
In this case both an algebra $A_{F}$ and a Hilbert space $H_{F}$ are finite
dimensional, and $D_{F}$ is just a matrix. One can go one step further and
consider products of spectral triples considered so far. The almost-
commutative geometry is a result of such a construction, where the first
spectral triple in the product is the canonical one for a manifold $M$, and
the other one is finite. In the case with $\dim M=4$, the Dirac operator for
the resulting triple is (pointwisely) $D_{M}\otimes 1+\gamma_{5}\otimes
D_{F}$, where $\gamma_{5}$ is the usual grading in the Clifford algebra
associated to the manifold $M$.

However, even for the product space with finite part being just the two points
set, this is not the most general Dirac operator one can consider. Indeed, the
operator

$\mathcal{D}=\begin{pmatrix}D_{1}&\gamma\Phi\\\
\gamma\Phi^{\ast}&D_{2}\end{pmatrix}$ (II.1)

with a field $\Phi$, which for our purposes is taken to be a constant, is an
example of another candidate sitarz2019 . Here $D_{1},D_{2}$ are two Dirac
operators for $M$, but considered with two different Riemannian metrics
$g_{1},g_{2}$. The operator $\gamma$ is a straightforward generalization of
$\gamma_{5}$: $\gamma^{\ast}=\gamma$, $\\{\gamma,\gamma^{a}\\}=0$ for all
anti-Hermitian $\gamma^{a}$ generating the Clifford algebra,
$\\{\gamma^{a},\gamma^{b}\\}=-2\delta^{ab}1$, but now $\gamma^{2}=\kappa=\pm
1$ (instead of requiring $\kappa=1$). These models are refered to as the
doubled geometries BS2021 .

We consider geometries of this type with the metric on each sheet chosen to be
of the form

$ds^{2}=\sum\limits_{j=0}^{3}a_{j}^{2}\left(dx^{j}\right)^{2},$ (II.2)

where $a_{j}$, for $j=0,1,2,3$, are constants. The spin connection $\omega$ is
identically zero since for the coframe $\\{\theta^{a}\\}$ we have
$d\theta^{a}=0$ for every $a=0,\ldots,3$, and the resulting Dirac operator
therefore reads,
$D=\sum\limits_{j=0}^{3}\frac{1}{a_{j}}\gamma^{j}\partial_{j}$. The
corresponding doubled geometry constructed out of these two sheets is
therefore described by a Dirac operator of the form

$\mathcal{D}=\sum\limits_{j=0}^{3}A_{j}\gamma^{j}\partial_{j}+\gamma F,$
(II.3)

where

$A_{j}=\begin{pmatrix}\frac{1}{a_{1,j}}&\\\
&\frac{1}{a_{2,j}}\end{pmatrix},\hskip 10.0ptF=\begin{pmatrix}&\Phi\\\
\Phi^{\ast}&\end{pmatrix}.$ (II.4)

The associated Laplace operator is hence given by

$\mathcal{D}^{2}=-\sum\limits_{j=0}^{3}A_{j}^{2}\partial_{j}^{2}+\sum\limits_{j=0}^{3}[F,A_{j}]\gamma\gamma^{j}\partial_{j}+\kappa
F^{2},$ (II.5)

and one can then easily read the decomposition of its symbol into the
homogeneous parts,
$\sigma_{\mathcal{D}^{2}}=\mathfrak{a}_{0}+\mathfrak{a}_{1}+\mathfrak{a}_{2}$.

Since our first goal is to determine the leading terms of the spectral action,

$\displaystyle\mathcal{S}(\mathcal{D})$
$\displaystyle=\Lambda^{4}\,\mathrm{Wres}(\mathcal{D}^{-4})+c\Lambda^{2}\,\mathrm{Wres}(\mathcal{D}^{-2})$
(II.6)
$\displaystyle=\int_{M}\int_{\|\xi\|=1}\left(\Lambda^{4}\,\mathrm{Tr}\,\mathrm{Tr}_{Cl}\,{\mathfrak{b}}_{0}^{2}+c\Lambda^{2}\,\mathrm{Tr}\,\mathrm{Tr}_{Cl}{\mathfrak{b}}_{2}\right),$

we have to find the symbol of the inverse of the Laplace operator,
$\sigma_{\mathcal{D}^{-2}}=\mathfrak{b}_{0}+\mathfrak{b}_{1}+\mathfrak{b}_{2}+...$,
what can be achieved by using the standard methods of pseudodifferential
calculus gilkey . (In the above equation $\mathrm{Tr}_{Cl}$ denotes the trace
performed over the Clifford algebra and $\mathrm{Tr}$ is the usual matrix
trace over two-by-two matrices.)

In our case we get

$\mathfrak{b}_{0}=\left(\sum\limits_{j=0}^{3}A_{j}^{2}\xi_{j}^{2}\right)^{-1},\qquad\mathfrak{b}_{2}=\mathfrak{b}_{0}\mathfrak{a}_{1}\mathfrak{b}_{0}\mathfrak{a}_{1}\mathfrak{b}_{0}-\mathfrak{b}_{0}\mathfrak{a}_{0}\mathfrak{b}_{0},$
(II.7)

so that

$\mathrm{Tr}_{Cl}(\mathfrak{b}_{2})=-4\kappa\mathfrak{b}_{0}\left(\sum\limits_{j=0}^{3}[F,A_{j}]\mathfrak{b}_{0}[F,A_{j}]\xi_{j}^{2}+F^{2}\right)\mathfrak{b}_{0}.$
(II.8)

The only nonzero elements of the matrix $\mathfrak{b}_{0}$ are on its diagonal
and they are equal to

$(\mathfrak{b}_{0})^{i}_{\
i}=\frac{1}{\sum\limits_{j=0}^{3}A_{i,j}^{2}\xi_{j}^{2}},$ (II.9)

where $A_{i,j}\equiv(A_{j})^{i}_{\ i}=\frac{1}{a_{i,j}}$, and as a result of a
straightforward computation we get

$\displaystyle\mathrm{Tr}\left(-4\kappa\mathfrak{b}_{0}\sum\limits_{j=0}^{3}[F,A_{j}]\mathfrak{b}_{0}[F,A_{j}]\mathfrak{b}_{0}\xi_{j}^{2}\right)$
(II.10)
$\displaystyle=4\kappa|\Phi|^{2}\sum\limits_{j,k=0}^{3}\frac{(A_{2,j}-A_{1,j})^{2}(A_{1,k}^{2}+A_{2,k}^{2})}{\left(\sum\limits_{l=0}^{3}A_{1,l}^{2}\xi_{l}^{2}\right)^{2}\left(\sum\limits_{l=0}^{3}A_{2,l}^{2}\xi_{l}^{2}\right)^{2}}\xi_{j}^{2}\xi_{k}^{2}.$

The resulting spectral action is therefore of the form

$\mathcal{S}(\mathcal{D})\sim\int_{M}\left(\Lambda_{e}^{2}\mathcal{S}_{\Lambda_{e}}+\alpha\widehat{V}(g_{1},g_{2})\right),$
(II.11)

with

$\mathcal{S}_{\Lambda_{e}}=\int_{\|\xi\|=1}\left\\{\left(\sum\limits_{j=0}^{3}A_{1,j}^{2}\xi_{j}^{2}\right)^{-2}+\left(\sum\limits_{j=0}^{3}A_{2,j}^{2}\xi_{j}^{2}\right)^{-2}\right\\}$
(II.12)

and

$\widehat{V}(g_{1},g_{2})=\sum\limits_{j,k=0}^{3}(A_{2,j}-A_{1,j})^{2}(A_{1,k}^{2}+A_{2,k}^{2})\int_{\|\xi\|=1}\frac{\xi_{j}^{2}\xi_{k}^{2}}{\left(\sum\limits_{l=0}^{3}A_{1,l}^{2}\xi_{l}^{2}\right)^{2}\left(\sum\limits_{l=0}^{3}A_{2,l}^{2}\xi_{l}^{2}\right)^{2}}$
(II.13)

where we have already introduced effective parametrization,

$\Lambda_{e}^{2}=\frac{12}{c}(\Lambda^{2}-c\kappa|\Phi|^{2}),\qquad\alpha=12|\Phi|^{2}\kappa,$
(II.14)

and ommited the irrelevant global multiplicative constant. Therefore, the
problem of finding the potential term describing the interaction between the
two diagonal metrics reduces to compute linear combination of the integrals of
the form

$\int_{\|\xi\|=1}\frac{\xi_{j}^{2}\xi_{k}^{2}}{\left(\sum\limits_{l=0}^{3}A_{1,l}^{2}\xi_{l}^{2}\right)^{2}\left(\sum\limits_{l=0}^{3}A_{2,l}^{2}\xi_{l}^{2}\right)^{2}}.$
(II.15)

Moreover, from the Eqn. (II.13) it immediately follows that
$\widehat{V}(g_{1},g_{2})=\widehat{V}(g_{2},g_{1})$, that is, $\widehat{V}$ is
symmetric under the interchange $g_{1}\leftrightarrow g_{2}$. We further
conjecture that the potential term can be written as

$\widehat{V}(g_{1},g_{2})=2\pi^{2}\mathbb{V}\left(\sqrt{g_{2}^{-1}g_{1}}\right)\sqrt{\det\,g_{2}}$
(II.16)

for some function $\mathbb{V}$. By the symmetry of $\widehat{V}$, to prove
this claim it is enough to show that the function
$\mathbb{V}^{\prime}(g_{1},g_{2}):=\frac{\widehat{V}(g_{1},g_{2})}{2\pi^{2}\sqrt{\det
g_{2}}}$ depends only on the eigenvalues of $\sqrt{g_{2}^{-1}g_{1}}$. We
illustrate this hypothesis on a simple nontrivial example - the Hopf model -
discussed in the forthcoming section.

##  III The Hopf model

We consider here models with diagonal metrics given by $g_{00}=g_{11}=b^{2}$
and $g_{22}=g_{33}=a^{2}$, for which then have

$(\mathfrak{b}_{0})^{1}_{\
1}=\frac{a_{1}^{2}b_{1}^{2}}{a_{1}^{2}(\xi_{0}^{2}+\xi_{1}^{2})+b_{1}(\xi_{2}^{2}+\xi_{3}^{2})},\qquad(\mathfrak{b}_{0})^{2}_{\
2}=\frac{a_{2}^{2}b_{2}^{2}}{a_{2}^{2}(\xi_{0}^{2}+\xi_{1}^{2})+b_{2}(\xi_{2}^{2}+\xi_{3}^{2})}$
(III.1)

and

$\displaystyle\widehat{V}(g_{1},g_{2})$
$\displaystyle=\frac{(b_{1}-b_{2})^{2}}{b_{1}^{2}b_{1}^{2}}\int_{\|\xi\|=1}(\xi_{0}^{2}+\xi_{1}^{2})\det(\mathfrak{b}_{0})\mathrm{Tr}(\mathfrak{b}_{0})$
(III.2)
$\displaystyle+\frac{(a_{1}-a_{2})^{2}}{a_{1}^{2}a_{1}^{2}}\int_{\|\xi\|=1}(\xi_{2}^{2}+\xi_{3}^{2})\det(\mathfrak{b}_{0})\mathrm{Tr}(\mathfrak{b}_{0}).$

In order to parametrize the three-sphere $\|\xi\|=1$ we use here the following
Hopf-like coordinates which resemble the symmetry of the system:

$\xi_{0}=\cos\theta\cos\varphi,\quad\xi_{1}=\cos\theta\sin\varphi,\quad\xi_{2}=\sin\theta\cos\psi,\quad\xi_{3}=\sin\theta\sin\psi.$
(III.3)

The angle $\theta$ is taken from $\left[0,\frac{\pi}{2}\right]$, while
$0\leq\psi\leq 2\pi$, and the surface element in these coordinates is then
given by $dS=\cos\theta\sin\theta d\theta d\varphi d\psi$.

As a result, we get

$\displaystyle\det(\mathfrak{b}_{0})\mathrm{Tr}(\mathfrak{b}_{0})$ (III.4)
$\displaystyle=$
$\displaystyle\frac{a_{1}^{2}a_{2}^{2}b_{1}^{2}b_{2}^{2}\left[a_{1}^{2}b_{1}^{2}(a_{2}^{2}\cos^{2}\theta+b_{2}^{2}\sin^{2}\theta)+a_{2}^{2}b_{2}^{2}(a_{1}^{2}\cos^{2}\theta+b_{1}^{2}\sin^{2}\theta)\right]}{\left[a_{1}^{2}\cos^{2}\varphi(a_{1}^{2}\cos^{2}\theta+b_{1}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}\left[a_{2}^{2}\cos^{2}\varphi(a_{2}^{2}\cos^{2}\theta+b_{2}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}}.$

Let us introduce the following notation

$I_{\mu,c}=\int_{\|\xi\|=1}\frac{\xi_{\mu}^{2}\cos^{2}\theta\,dS}{\left[a_{1}^{2}\cos^{2}\varphi(a_{1}^{2}\cos^{2}\theta+b_{1}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}\left[a_{2}^{2}\cos^{2}\varphi(a_{2}^{2}\cos^{2}\theta+b_{2}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}},$
(III.5)

and

$I_{\mu,s}=\int_{\|\xi\|=1}\frac{\xi_{\mu}^{2}\sin^{2}\theta\,dS}{\left[a_{1}^{2}\cos^{2}\varphi(a_{1}^{2}\cos^{2}\theta+b_{1}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}\left[a_{2}^{2}\cos^{2}\varphi(a_{2}^{2}\cos^{2}\theta+b_{2}^{2}\sin^{2}\theta)\sin^{2}\varphi\right]^{2}},$
(III.6)

and notice that

$I_{0,c}=I_{1,c}=I_{2,c}=I_{3,c},\qquad I_{0,s}=I_{1,s},\qquad
I_{2,s}=I_{3,s},$ (III.7)

so that

$\displaystyle\widehat{V}(g_{1},g_{2})$
$\displaystyle=2a_{1}^{2}a_{2}^{2}b_{1}^{2}b_{2}^{2}\left\\{a_{1}^{2}a_{2}^{2}(b_{1}^{2}+b_{2}^{2})\left[\frac{(b_{1}-b_{2})^{2}}{b_{1}^{2}b_{2}^{2}}I_{0,c}+\frac{(a_{1}-a_{2})^{2}}{a_{1}^{2}a_{2}^{2}}I_{0,s}\right]\right.$
(III.8)
$\displaystyle\left.+b_{1}^{2}b_{2}^{2}(a_{1}^{2}+a_{2}^{2})\left[\frac{(b_{1}-b_{2})^{2}}{b_{1}^{2}b_{2}^{2}}I_{0,s}+\frac{(a_{1}-a_{2})^{2}}{a_{1}^{2}a_{2}^{2}}I_{2,s}\right]\right\\}.$

In order to find the final form of the potential it remains to compute the
integrals $I_{0,c}$, $I_{0,s}$ and $I_{2,s}$. The result reads,

$\widehat{V}(g_{1},g_{2})=\frac{2\pi^{2}}{(a_{2}b_{1}-a_{1}b_{2})(a_{2}b_{1}+a_{1}b_{2})^{2}}\left(F(a_{1},a_{2},b_{1},b_{2})+G(a_{1},a_{2},b_{1},b_{2})\right),$
(III.9)

where

$F(a_{1},a_{2},b_{1},b_{2})=4a_{1}^{2}a_{2}^{2}b_{1}^{2}b_{2}^{2}(a_{1}-a_{2})(b_{1}-b_{2})\log\left(\frac{a_{1}b_{2}}{a_{2}b_{1}}\right),$
(III.10)

and

$\displaystyle
G(a_{1},a_{2},b_{1},b_{2})=(a_{2}^{2}b_{1}^{2}-a_{1}^{2}b_{2}^{2})\left[a_{1}^{2}b_{1}^{2}a_{2}(b_{1}-2b_{2})\right.$
$\displaystyle\left.+a_{2}^{2}b_{2}^{2}a_{1}(b_{2}-2b_{1})\right.$ (III.11)
$\displaystyle+\left.a_{1}^{3}b_{1}^{2}b_{2}+a_{2}^{3}b_{2}^{2}b_{1}\right].$

We observe that for $a_{1}=a_{2}$ the potential reduces to
$\widehat{V}(g_{1},g_{2})=2\pi^{2}a_{1}^{2}(b_{1}-b_{2})^{2}$, while for
$b_{1}=b_{2}$ we have
$\widehat{V}(g_{1},g_{2})=2\pi^{2}(a_{1}-a_{2})^{2}b_{1}^{2}$.

Since the logarithm vanishes if and only if
$\frac{a_{1}}{a_{2}}=\frac{b_{1}}{b_{2}}$ it would be, in principle,
interesting to consider the limit of $\widehat{V}(g_{1},g_{2})$ when $b_{1}$
tends to $b_{2}\frac{a_{1}}{a_{2}}$. The value of the function $V$ is
indeterminated in this case, but the limit may still exists. Indeed, as a
result we get

$\lim_{b_{1}\rightarrow
b_{2}\frac{a_{1}}{a_{2}}}\widehat{V}(g_{1},g_{2})=\frac{b_{2}^{2}}{a_{2}^{2}}(a_{1}-a_{2})^{2}(a_{1}^{2}+a_{2}^{2}).$
(III.12)

Introducing the new variable $x=\frac{b_{1}}{b_{2}}$ and
$y=\frac{a_{1}}{a_{2}}$ we can write

$\widehat{V}(g_{1},g_{2})=2\pi^{2}\mathbb{V}\left(\sqrt{g_{2}^{-1}g_{1}}\right)\sqrt{\det\,g_{2}},$
(III.13)

where

$\mathbb{V}\left(\sqrt{g_{2}^{-1}g_{1}}\right)=\frac{4x^{2}y^{2}(x-1)(y-1)}{(x-y)(x+y)^{2}}\log\left(\frac{y}{x}\right)+x^{2}y^{2}+1-2xy\frac{xy+1}{x+y}.$
(III.14)

We observe that

$\mathbb{V}\left(\sqrt{g_{2}^{-1}g_{1}}\right)\sqrt{\det\,g_{2}}=\mathbb{V}\left(\sqrt{g_{1}^{-1}g_{2}}\right)\sqrt{\det\,g_{1}},$
(III.15)

what illustrates the hypothesis.

##  IV Comment on generic metrics

The so far examined examples of doubled geometries suggest that these models
can be thought of as certain modifications of bimetric theories as the
potential term possesses features characteristic to this type of modified
gravity theories. Despite the fact that series of non-trivial examples are
already analysed, the derivation of the action in the generic case is still an
open problem. In the approach we are using the main chalenge is related with
the computation of certain integrals of rational functions defined over higher
spheres:

$I=\int_{\|\xi\|=1}\frac{d^{4}\xi}{A_{\mu\nu}\xi^{\mu}\xi^{\nu}},$ (IV.1)

with smooth $A_{\mu\nu}$, which can be further written as
$A_{\mu\nu}=\Omega(\delta_{\mu\nu}+\epsilon_{\mu\nu})$ with
$\Omega\in\mathbb{R}$ and $\epsilon_{\mu\nu}$ being symmetric and traceless.
In the formula above $\xi^{\alpha}$ is the $\alpha$th coordinate of vector
$\xi$.

We make here some comments on the analysis of doubled geometry models in case
when the tensors $A_{\mu\nu}$ does not differ sufficiently from the diagonal
ones. This is not identical to the situation where the metrics are small
perturbation of the Euclidean ones - for the discussion of the latter we refer
to BoSi22 .

Assuming that $\max\limits_{\mu,\nu}\|\epsilon_{\mu\nu}\|_{\infty}$ is
sufficiently small we can expand in these parameters and write

$I=\frac{1}{\Omega}\sum\limits_{m\geq
0}(-2)^{m}\epsilon_{\alpha_{1}\beta_{1}}\ldots\epsilon_{\alpha_{m}\beta_{m}}I^{\alpha_{1}\beta_{1}\ldots\alpha_{m}\beta_{m}}_{4,m},$
(IV.2)

where

$I^{\alpha_{1}\beta_{1}\ldots\alpha_{m}\beta_{m}}_{n,m}=\int_{\mathbb{S}^{n-1}}d^{n}\xi\xi^{\alpha_{1}}\xi^{\beta_{1}}\ldots\xi^{\alpha_{m}}\xi^{\beta_{m}}$
(IV.3)

are polynomial integrals over higher spheres which can be evaluated
generalizing the methods from othmani  \- see also BoSi22  for further
discussion.

Denoting

$\gamma_{j}=\begin{cases}\alpha_{k},\quad j=2k-1,\\\ \beta_{k},\quad
j=2k\end{cases}$ (IV.4)

we define
$I^{\gamma_{1}\ldots\gamma_{2m}}=I^{\alpha_{1}\beta_{1}\ldots\alpha_{m}\beta_{m}}_{4,m}$,
and let $\Delta^{\gamma_{1}...\gamma_{2m}}$ be the sum of product of deltas in
$I^{\gamma_{1}...\gamma_{2m}}$, i.e.

$I^{\gamma_{1}...\gamma_{2m}}=c_{m}\Delta^{\gamma_{1}...\gamma_{2m}}$ (IV.5)

with some real number $c_{m}$, and
$\Delta^{...}=\sum\delta^{..}...\delta^{..}$ Since $\epsilon$ is traceless not
all terms in

$\epsilon_{\gamma_{1}\gamma_{2}}...\epsilon_{\gamma_{2m-1}\gamma_{2m}}I^{\gamma_{1}...\gamma_{2m}}$
(IV.6)

are nonzero. Let $\mathcal{N}_{2m}$ be the number of nonzero terms, and
consider the following problem. Suppose the numbers $1,...,2m$ are given, and
we would like to use them to fill in an $1\times 2m$ array $T$, with a given
subdivision into $1\times 2$ subarrays $T=T_{1}|T_{2}|...|T_{m}$, as follows:

  * •

In the first entry of $T_{1}$ we put $1$,

  * •

For every $j=1,...,m$, we have $T_{j}=\begin{bmatrix}a_{j}|b_{j}\end{bmatrix}$
with $a_{j}<b_{j}$,

  * •

For every $j=1,...,m$, $a_{j}<a_{j+1}$.

Then $\mathcal{N}_{2m}$ is a number of such fillings for which there is no $j$
such that $T_{j}$ is of the form $\begin{bmatrix}2l-1|2l\end{bmatrix}$, for
some $l=1,...,m$.

Since, by symmetry of $\epsilon$, any nonzero term in
$\epsilon_{\gamma_{1}\gamma_{2}}...\epsilon_{\gamma_{2m-1}\gamma_{2m}}I^{\gamma_{1}...\gamma_{2m}}$
produces $\mathrm{tr}(\epsilon^{m})$, we have

$\epsilon_{\gamma_{1}\gamma_{2}}...\epsilon_{\gamma_{2m-1}\gamma_{2m}}I^{\gamma_{1}...\gamma_{2m}}=\mathcal{N}_{2m}c_{m}\mathrm{tr}(\epsilon^{m}),$
(IV.7)

and the problem reduces to finding coefficients $c_{m}$. Since
$\mathrm{area}(\mathbb{S}^{3})=2\pi^{2}$ we get $c_{1}=\frac{\pi^{2}}{2}$.
Moreover, by using the generalization of (othmani, , Prop. 2) (see also
(BoSi22, , Prop. A.2) ) one can easily find the recursive formula for $c_{m}$:

$c_{m}=\frac{c_{m-1}}{4+2(m-1)},$ (IV.8)

and its solution reads

$c_{m}=\frac{4\pi^{2}}{(2m+2)!!}.$ (IV.9)

As a result,

$I=\frac{1}{\Omega}\left[2\pi^{2}+\frac{2\pi^{2}}{3}\mathrm{tr}(\epsilon^{2})+4\pi^{2}\sum_{m\geq
3}\frac{(-2)^{m}}{(2m+2)!!}\mathcal{N}_{2m}\mathrm{tr}(\epsilon^{m})\right].$
(IV.10)

In order to apply this result to a specific term of the action one has to
first solve the combinatorial problem of finding the coefficients
$\mathcal{N}_{2m}$, up to required order in $m$. We postpone for the future
research the problem of determining set of metrics for which the rate of
convergence of the above series is satisfactory for all the terms that appear
in the action functional.

##  V Conclusions and outlook

The discussed doubled geometry model is an interesting possibility of going
beyond the General Relativity. The explicit functional form of its action is
derivable in the same way as the Hilbert-Einstein’s one but with the use of a
different geometry instead of the classical manifold. Here we extended the
existing family of known examples for which the features characteristic to
bimetric gravity models are present. We also made further steps towards the
analysis of models that are beyond the class of such whose action is
analytically computable. We remark that yet another approach based on a
different type of noncommutativity can produce bimetric type of models Ce18 .
It will be interesting to find some deeper relations between these two
formulations - we postpone this for a future research.

Acknowledgments
The author thanks A. Sitarz for helpful discussions and comments. The author
acknowledges the support from the National Science Centre, Poland, Grant No.
2018/31/N/ST2/00701.

## References

  * (1) A. Connes, Noncommutative Geometry, Academic Press, 1994.
  * (2) A. Connes, On the spectral characterization of manifolds, J. Noncommut. Geom. 7, 1–82 (2013).
  * (3) A. Connes, Noncommutative geometry and reality, J. Math. Phys. (N.Y.) 36, 6194 (1995).
  * (4) A. Connes, Gravity coupled with matter and foundation of non-commutative geometry, Commun. Math. Phys. 182, 155 (1996).
  * (5) A. H. Chamseddine, A. Connes and M. Marcolli, Gravity and the standard model with neutrino mixing, Adv. Theor. Math. Phys. 11, 991 (2007).
  * (6) F. Lizzi, Noncommutative Geometry and Particle Physics, Proceedings of Corfu Summer Institute 2017 “Schools and Workshops on Elementary Particle Physics and Gravity” - PoS(CORFU2017) (2018).
  * (7) S. F. Hassan and R. A. Rosen, Bimetric gravity from ghost-free massive gravity, J. High Energ. Phys. 02, 126 (2012).
  * (8) Y. Akrami, T. S. Koivisto, D. F. Mota and M. Sandstad, Bimetric gravity doubly coupled to matter: theory and cosmological implications, J. Cosmol. Astropart. Phys. 10 (2013) 046.
  * (9) W. D. van Suijlekom, Noncommutative Geometry and Particle Physics, Springer Netherlands, Dordrecht 2015.
  * (10) A. Connes and J. Lott, Particle models and noncommutative geometry, Nucl. Phys. B, Proc. Suppl. 18, 29 (1991).
  * (11) T. Krajewski, Classification of finite spectral triples, J. Geom. Phys. 28, 1-30 (1998).
  * (12) M. Paschke and A. Sitarz, Discrete spectral triples and their symmetries, J. Math. Phys. 39, 6191 (1998).
  * (13) A. Bochniak and A. Sitarz, Stability of Friedmann-Lemaître-Robertson-Walker solutions in doubled geometries, Phys. Rev. D 103 (2021) 044041.
  * (14) A. Sitarz, On almost commutative Friedmann-Lemaître-Robertson-Walker geometries, Class. Quantum Grav. 36 (2019) 195007.
  * (15) P. Gilkey, Invariance theory, the heat equation, and the Atiyah-Singer index theorem, Mathematics Lecture Series 11. Publish or Perish 1984.
  * (16) S. F. Hassan and R. A. Rosen,Bimetric gravity from ghost-free massive gravity, J. High Energ. Phys. 2012, 126 (2012)
  * (17) Y. Akrami, T. S. Koivisto, D. F. Mota and M. Sadstad, Bimetric gravity doubly coupled to matter: theory and cosmological implications, Journal of Cosmology and Astroparticle Physics, 10 (2013) 046.
  * (18) Y. Othmani, Polynomial integration over the unit sphere, Applied Mathematics Letters 24(7)(2011) 1260–1264.
  * (19) A. Bochniak and A. Sitarz, Spectral interaction between universes, arXiv:2201.03839.
  * (20) M. de Cesare, M. Sakellariadou, and P. Vitale, Noncommutative gravity with self-dual variables, Classical Quantum Gravity 35, 215009 (2018).