# Causal games of work extraction with indefinite causal order

Gianluca Francica<EMAIL_ADDRESS>Dipartimento di Fisica e
Astronomia “G. Galilei”, Università degli Studi di Padova, via Marzolo 8,
35131 Padova, Italy

###### Abstract

An indefinite causal order, where the causes of events are not necessarily in
past events, is predicted by the process matrix framework. A fundamental
question is how these non-separable causal structures can be related to the
thermodynamic phenomena. Here, we approach this problem by considering the
existence of two cooperating local Maxwell’s demons which try to exploit the
presence of global correlations and indefinite causal order to optimize the
extraction of work. Thus, we prove that it is possible to have a larger
probability to lower the local energy to zero if causal inequalities are
violated, and that can be extracted more average work with respect to a
definite causal order. However, for non-interacting parties, for the system
considered the work extractable cannot be larger than the definite causal
order bound.

##  I Introduction

In the common view of the nature, the events have a definite causal order,
i.e., the causes of events are in past events. In contrast, it is known that
the process matrix formalism of Oreshkov, Costa and Brukner Oreshkov12  allows
causal structures compatible with local quantum mechanics for two parties
Alice and Bob, which are causally non-separable, i.e., neither Alice comes
before Bob nor Bob comes before Alice, nor a mixture thereof. This indefinite
causal order can be exploited in some “causal games”, e.g., allows us to find
a strategy violating the so-called causal inequalities Oreshkov12 ;
Branciard16 , even if there also are causally non-separable processes which
admit a causal model Feix16 , i.e., they do not violate causal inequalities (a
popular example is the quantum switch model Chiribella13 ). Concerning the
question if there is causation in fundamental physics, some relevant insights
are recently given in Ref. Adlam22 . A further fundamental question, which
here we take in exam, can be if indefinite causal order can be compatible with
the second law of thermodynamics. Furthermore, concerning the applications, we
note that indefinite causal order can provide advantages in certain
communication feix15 ; guerin16 ; ebler18  and computing tasks chiribella12 ;
araujo14 , and recent experiments have been performed procopio15 ; rubino17 ;
goswami18 ; goswami20 ; rubino21 . Recently, a certain attention has also been
paid to understand the role of indefinite causal order to achieve
thermodynamic tasks Nie20 ; Cao21 ; Felce21 ; Felce20 ; Guha20 ; Simonov22 ;
Dieguez22 ; Guha22 ; Capela22 . We note that the daemonic ergotropy Francica17
, where the information gained by performing measurements on a part of the
system can be communicated to locally extract work from a different part, can
be naturally related to a causal structure. Basically, in the daemonic
ergotropy extraction Alice and Bob share a correlated state, if Alice performs
measurements on her part and communicates to Bob the outcomes, Bob can perform
an optimal unitary cycle to extract the maximum work from his part. Of course,
the measurements performed by Alice are in the past of Bob, thus the protocol
exhibits a definite causal order. Here, we adopt the same point of view, with
the addition that the roles of Alice and Bob can be also reversed, in order to
get a causally non-separable structure. As we will show, this allows us to get
a work extraction where indefinite causal order plays some major role.

##  II Preliminaries of thermodynamics

Let us introduce some concepts of work extraction in finite quantum systems.
We are interested in extracting work from a quantum system in an initial state
$\rho_{0}$ by using a cyclic unitary transformation, i.e., by considering a
unitary time evolution operator $U$ generated by a time-dependent Hamiltonian
which at the final time is equal to the initial Hamiltonian $H$. We aim to
maximize the average work extracted from the system, or equivalently to
minimize the final average energy $\langle
E\rangle=\text{Tr}\left\\{H\rho\right\\}$, where $\rho=U\rho_{0}U^{\dagger}$
is the final state, over the set of all the unitary operators $U$, since in
general the average work extracted is equal to minus the change of average
energy if the system is thermally isolated. As shown in Ref. allahverdyan04 ,
the final average energy $\langle E\rangle$ is minimal if the final state is a
passive state, i.e.,
$\rho=P_{\rho_{0}}\equiv\sum_{k}r_{k}\left|\epsilon_{k}\right\rangle\left\langle\epsilon_{k}\right|$,
where we have ordered the labels of eigenstates of $H$ and of $\rho_{0}$ such
that
$H=\sum_{k}\epsilon_{k}\left|\epsilon_{k}\right\rangle\left\langle\epsilon_{k}\right|$,
with $\epsilon_{k}\leq\epsilon_{k+1}$, and
$\rho_{0}=\sum_{k}r_{k}\left|r_{k}\right\rangle\left\langle r_{k}\right|$,
with $r_{k}\geq r_{k+1}$. The resulting average work is known as ergotropy. In
particular, for an initial pure state, the passive state will be the ground-
state $\left|\epsilon_{1}\right\rangle$. Concerning the local extraction of
work, if the system consists of two parties, Alice and Bob, which do not
interact with each other, Alice can perform measurements on her part,
described by a set of orthogonal projectors $\Pi^{A}_{x}$ of rank one, and can
communicate the outcomes $x$ to Bob. Thus, the state of Bob collapses into
$\rho_{B|x}=\text{Tr}_{A}\left\\{\Pi^{A}_{x}\otimes
I^{B}\rho_{0}\right\\}/p(x)$ with probability
$p(x)=\text{Tr}\left\\{\Pi^{A}_{x}\otimes I^{B}\rho_{0}\right\\}$, where
$I^{B}$ is the identity matrix on the Bob’s Hilbert space, and Bob can perform
a local cyclic unitary transformation $U_{x}$ conditioned on the outcome $x$
of the measurement. We aim to maximize the average work extracted from the
Bob’s part, or equivalently to minimize the final average energy of Bob
$\langle
E_{B}\rangle=\sum_{x}p(x)\text{Tr}\left\\{H_{B}U_{x}\rho_{B|x}U_{x}^{\dagger}\right\\}$
over the set of all the unitary operators $U_{x}$ for given projectors
$\Pi^{A}_{x}$, where $H_{B}$ is the Hamiltonian of Bob. The resulting average
work extracted by Bob, known as daemonic ergotropy, is minus the change of the
average energy of Bob and has been investigated in terms of the correlations
between the two parties in Ref. Francica17 . By considering that the final
average energy of Alice is $\langle
E_{A}\rangle=\sum_{x}\text{Tr}\left\\{\rho_{0}\Pi^{A}_{x}H_{A}\Pi^{A}_{x}\otimes
I^{B}\right\\}$, where $H_{A}$ is the Hamiltonian of Alice, the final average
energy of the total system is $\langle E\rangle=\langle E_{A}\rangle+\langle
E_{B}\rangle$. Typically, the average energy of Alice can change because of
the measurements, however for the special case where the Hilbert spaces of
Alice and Bob have the same dimension $d$ and $\rho_{0}$ is a maximally
entangled pure state, the reduced states of Alice and Bob are equal to the
completely mixed state, so that the average energy of Alice does not change
for any measurements and the daemonic ergotropy, i.e., the average work
extracted $\langle w\rangle$ is minus the change of average energy of the
total system. In particular, for any projective measurements $\Pi^{A}_{x}$, we
get $\langle w\rangle=\sum_{i}\epsilon^{B}_{i}/d-\epsilon^{B}_{1}$, where
$\epsilon^{B}_{i}$ are the eigenvalues of $H_{B}$ and $\epsilon^{B}_{1}$ is
the lowest one, because Bob can perform local cyclic unitary transformations
$U_{x}$ to always achieve the state with lowest energy. Of course, in this
scheme of work extraction, we have a definite causal order. In the next
section, with the aim to consider causally non-separable structures, we will
modify this scheme. We will be interested in the average work extracted and so
in the final average energy $\langle E\rangle$ of the total system, which,
given the final state $\rho$, can be calculated as $\langle
E\rangle=\text{Tr}\left\\{\rho H\right\\}$. Moreover, we will be interested in
a certain subsystem $X$, having final average energy $\langle
E_{X}\rangle=\text{Tr}\left\\{\rho_{X}H_{X}\right\\}$, where $\rho_{X}$ and
$H_{X}$ are the final reduced state and the Hamiltonian of the subsystem,
respectively. In particular, we note that in general the final energy of the
subsystem $X$ has the probability distribution
$p(E_{X})=\sum_{i}\left\langle\epsilon^{X}_{i}\right|\rho_{X}\left|\epsilon^{X}_{i}\right\rangle\delta(E_{X}-\epsilon^{X}_{i})$,
where $\left|\epsilon^{X}_{i}\right\rangle$ is the eigenstate of $H_{X}$ with
eigenvalue $\epsilon^{X}_{i}$, so that $\langle E_{X}\rangle=\int
E_{X}p(E_{X})dE_{X}$.

##  III Causal games of work extraction

To link the daemonic ergotropy scheme to an indefinite causal order structure,
we assume that Alice and Bob share two correlated states. In particular both
Alice and Bob have a square qubit and a circle qubit. The square qubit of
Alice and the circle one of Bob (red qubits) are correlated, in particular are
in a singlet state, and the circle qubit of Alice and the square one of Bob
(blue qubits) are also in a singlet state. Alice and Bob can perform local
operations and can communicate each other. In detail, Alice and Bob perform
measurements on the square qubits, and unitary operations on the circle ones.
A schematic illustration of the system is given in Fig. 1.

Figure 1: The system is made of four parties, represented by two squares and
two circles. The parties in $A$ and $B$ belong to Alice and Bob, respectively,
and the parties connected by a line are initially correlated. Local
measurements and unitary operations are performed on the squares and the
circles, respectively.

Since Alice and Bob can communicate, due to the presence of global
correlations, in principle they can perform optimized unitary operations
depending on the measurements to achieve a precise task. Here, we consider the
possibility to have causally non-separable processes, and a thermodynamic task
that is the work extraction from the qubits system, i.e., from the two singlet
states. What is the role of indefinite causal order in achieving this task?
Due to the violation of causal inequalities, in particular, the probability to
guess the neighbor values of the measurement is larger in the causally non-
separable case Branciard16 , one can expect to extract more work than in the
causally separable case. However, the answer is not that simple, let us show
why. We start to define the Hamiltonian of the system. The qubits do not
interact each other, and a qubit has the Hamiltonian
$H_{q}=\epsilon(\sigma_{x}+1)/2$, with $\epsilon>0$, where $\sigma_{\alpha}$,
with $\alpha=x,y,z$, are the Pauli matrices. Alice and Bob will perform
measurements of $\sigma_{z}$ on their square qubits (the main results do not
change by changing basis). If the system is in the initial state, Alice gets
the bit $x$ and the singlet state collapses in
$\left|x\right\rangle\otimes\left|\bar{x}\right\rangle$ with probability
$p(x)=1/2$, where $x=0,1$, $\left|x\right\rangle$ is eigenstate of
$\sigma_{z}$ with eigenvalue $2x-1$ and $\bar{x}=0$ if $x=1$, $\bar{x}=1$ if
$x=0$. Similarly, Bob gets the bit $y$ and the singlet state collapses in
$\left|y\right\rangle\otimes\left|\bar{y}\right\rangle$ with probability
$p(y)=1/2$. We note that one measurement does not change the average energy of
a singlet state which is $\epsilon$ for any value of the bit obtained. This
energy can be lowered due to a local cyclic unitary transformation. We start
to take in exam a definite causal order where Alice precedes Bob, i.e.,
$A\prec B$, so that, at first, Alice performs a measurement on her square red
qubit, and communicate to Bob the bit $x$, and performs a local cyclic unitary
transformation $U_{A}$ on her circle blue qubit. Later, Bob uses the bit $x$
of information received to perform a local cyclic unitary transformation
$U_{x}$ on his circle red qubit, so that the initial red singlet state becomes
$\left|x\right\rangle\otimes U_{x}\left|\bar{x}\right\rangle$ with probability
$p(x)=1/2$, thus with the aim to lower the energy, Bob performs the unitary
$U_{x}$ such that $U_{x}\left|\bar{x}\right\rangle=\left|-\right\rangle$,
where $\left|\pm\right\rangle$ is the eigenstate of $\sigma_{x}$ with
eigenvalue $\pm 1$, and the final average energy of the two qubits will be
$\epsilon/2$ with equal probability, so that its average is $\langle
E_{red}\rangle=\epsilon/2$. In detail, the local unitary $U_{x}$ is such that
$U_{0}\left|1\right\rangle=\left|-\right\rangle$,
$U_{0}\left|0\right\rangle=\left|+\right\rangle$,
$U_{1}\left|0\right\rangle=\left|-\right\rangle$ and
$U_{1}\left|1\right\rangle=\left|+\right\rangle$. The blue qubits are in the
state $(U_{A}\left|0\right\rangle\otimes\left|1\right\rangle-
U_{A}\left|1\right\rangle\otimes\left|0\right\rangle)/\sqrt{2}$ with an
average energy equal to $\epsilon$ for any $U_{A}$. If Bob performs a
measurement on his square blue qubit, we get the state
$U_{A}\left|\bar{y}\right\rangle\otimes\left|y\right\rangle$ with probability
$p(y)=1/2$, but in average the energy remains $\langle
E_{blue}\rangle=\epsilon$. Thus, in this case the final average energy of the
four qubits is equal to $\langle E\rangle=\langle E_{red}\rangle+\langle
E_{blue}\rangle=\epsilon/2+\epsilon=3\epsilon/2$, which is the lowest value
gettable for a definite causal order. Since the measurements do not change the
average energy, the average work extracted is $\langle
w\rangle=2\epsilon-\langle E\rangle$ and gets the maximum value $\langle
w\rangle=\epsilon/2$, which is also equal to the daemonic ergotropy of the red
qubits. Of course, if Alice and Bob do not communicate but only use local
operations, in average they cannot lower the energy of the singlet states, so
that the average work is zero, $\langle w\rangle=0$, in agreement with the
second law of thermodynamics. Thus, the gain achieved is related to the use of
the information acquired.

In contrast, in the presence of indefinite causal order Alice and Bob are
given the bit inputs $x$ and $y$, with probability $p(x,y)$, and return the
bit outputs $a$ and $b$, respectively. By taking in account Alice, for each
input $x$ and output $a$, we associate an operation described by a completely
positive map
$\mathcal{M}^{A_{I}A_{O}}_{a|x}:\mathcal{L}(\mathcal{H}^{A_{I}})\to\mathcal{L}(\mathcal{H}^{A_{O}})$,
where $\mathcal{L}(\mathcal{H}^{X})$ is the space of linear operators over the
Hilbert space $\mathcal{H}^{X}$ of dimension $d_{X}=2$. We note that all the
maps must sum up to a trace-preserving map. Using the Choi-Jamiolłkowski
isomorphism Choi75 ; Jamiolkowski72 , we represent the map
$\mathcal{M}^{A_{I}A_{O}}_{a|x}$ as the operator
$M^{A_{I}A_{O}}_{a|x}=[I^{A_{I}}\otimes\mathcal{M}^{A_{I}A_{O}}_{a|x}(\left|\varphi^{+}\right\rangle\left\langle\varphi^{+}\right|)]^{T}\in\mathcal{L}(\mathcal{H}^{A_{I}}\otimes\mathcal{H}^{A_{O}})$,
where $I^{X}$ is the identity matrix on $\mathcal{H}^{X}$ and
$\left|\varphi^{+}\right\rangle=\sum_{i}\left|ii\right\rangle$. The operators
$M^{A_{I}A_{O}}_{a|x}$ are such that $M^{A_{I}A_{O}}_{a|x}\geq 0$ for each $a$
and
$\text{Tr}_{A_{O}}\left\\{\sum_{a}M^{A_{I}A_{O}}_{a|x}\right\\}=I^{A_{I}}$.
Similarly, for Bob we get the operators
$M^{B_{I}B_{O}}_{b|y}\in\mathcal{L}(\mathcal{H}^{B_{I}}\otimes\mathcal{H}^{B_{O}})$.
The joint conditional probability reads

$p(a,b|x,y)=\text{Tr}\left\\{(M^{A_{I}A_{O}}_{a|x}\otimes
M^{B_{I}B_{O}}_{b|y})W\right\\}\,,$ (1)

where
$W\in\mathcal{L}(\mathcal{H}^{A_{I}}\otimes\mathcal{H}^{A_{O}}\otimes\mathcal{H}^{B_{I}}\otimes\mathcal{H}^{B_{O}})$
is the so-called process matrix, which is an hermitian operator such that the
probabilities given by Eq. (1) are non-negative and normalized. In particular,
a process matrix $W$ needs to satisfy the conditions Araujo15

$\displaystyle W$ $\displaystyle\geq$ $\displaystyle 0\,,$ (2)
$\displaystyle\text{Tr}\left\\{W\right\\}$ $\displaystyle=$ $\displaystyle
d_{A_{O}}d_{B_{O}}\,,$ (3) ${}_{B_{I}B_{O}}W$ $\displaystyle=$
${}_{A_{O}B_{I}B_{O}}W\,,$ (4) ${}_{A_{I}A_{O}}W$ $\displaystyle=$
${}_{A_{I}A_{O}B_{O}}W\,,$ (5) $\displaystyle W$ $\displaystyle=$
${}_{B_{O}}W+{{}_{A_{O}}W}-{{}_{A_{O}B_{O}}W}\,,$ (6)

where we have defined the operation

${}_{X}W=\frac{I^{X}}{d_{X}}\otimes\text{Tr}_{X}\left\\{W\right\\}\,.$ (7)

If Bob cannot signal to Alice or Alice cannot signal to Bob we have the
process matrices $W^{A\prec B}=W^{A_{I}A_{O}B_{I}}\otimes I^{B_{O}}$ or
$W^{B\prec A}=W^{A_{I}B_{I}B_{O}}\otimes I^{A_{O}}$, respectively, where
$W^{A_{I}A_{O}B_{I}}\in\mathcal{L}(\mathcal{H}^{A_{I}}\otimes\mathcal{H}^{A_{O}}\otimes\mathcal{H}^{B_{I}})$
and
$W^{A_{I}B_{I}B_{O}}\in\mathcal{L}(\mathcal{H}^{A_{I}}\otimes\mathcal{H}^{B_{I}}\otimes\mathcal{H}^{B_{O}})$.
In detail, from Eqs. (2-6) we get the conditions for the matrix
$W^{A_{I}A_{O}B_{I}}$

$\displaystyle W^{A_{I}A_{O}B_{I}}$ $\displaystyle\geq$ $\displaystyle 0\,,$
(8) $\displaystyle\text{Tr}\left\\{W^{A_{I}A_{O}B_{I}}\right\\}$
$\displaystyle=$ $\displaystyle d_{A_{O}}\,,$ (9)
${}_{B_{I}}W^{A_{I}A_{O}B_{I}}$ $\displaystyle=$
${}_{A_{O}B_{I}}W^{A_{I}A_{O}B_{I}}\,,$ (10)

and similar conditions for the matrix $W^{A_{I}B_{I}B_{O}}$. Therefore, a
process is causally separable if the process matrix can be expressed in a
convex combination

$W_{sep}=qW^{A\prec B}+(1-q)W^{B\prec A}\,.$ (11)

We recall that causal non-separability can be inferred by using causal
inequalities Oreshkov12 ; Branciard16 ; Araujo15 ; Baumeler14 ; Oreshkov16 .
An example of causal inequality is a bound of the probability of success of
the “guess your neighbor’s input” game Branciard16 . For uniform input bits
$x$ and $y$, the probability of success is
$p_{succ}=1/4\sum_{x,y}p(a=y,b=x|x,y)$ and for a separable process
$p_{succ}\leq 1/2$, but it is known that there are causally non-separable
processes such that $p_{succ}>1/2$. Anyway, there also are non-separable
process, having a causal model, that do not violate causal inequalities Feix16
. One can expect that $p_{succ}$ will play a role in the work extraction. The
work extraction scheme can be generalized to the case of indefinite causal
order by requiring that Alice and Bob perform the local cyclic unitary
transformations $U_{a}$ and $U_{b}$, respectively, which depend on the bits
$a$ and $b$. It is easy to see that the final state is the mixture

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{x,y,a,b}p(x,y)p(a,b|x,y)\left|x\right\rangle\left\langle
x\right|\otimes
U_{b}\left|\bar{x}\right\rangle\left\langle\bar{x}\right|U_{b}^{\dagger}$ (12)
$\displaystyle\otimes
U_{a}\left|\bar{y}\right\rangle\left\langle\bar{y}\right|U_{a}^{\dagger}\otimes\left|y\right\rangle\left\langle
y\right|\,,$

where in our case $p(x,y)=1/4$. Of course the two initial non-correlated
singlets become correlated but are still separable. The final average energy
of the state $\rho$ is

$\langle
E\rangle=\epsilon+p_{1}\epsilon+2p_{2}\epsilon=2\epsilon-(p_{succ}-p_{2})\epsilon\,,$
(13)

where we have defined the probabilities $p_{1}$ and $p_{2}$ to wrong one and
two bits, respectively, which are
$p_{2}=1/4\sum_{x,y}p(a=\bar{y},b=\bar{x}|x,y)$ and $p_{1}=1-p_{succ}-p_{2}$.
To derive Eq. (13), we observe that the final state $\rho$ can be expressed as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{x,y}\big{(}p(a=y,b=x|x,y)P(\left|x\right\rangle\otimes\left|-\right\rangle\otimes\left|-\right\rangle\otimes\left|y\right\rangle)$
(14)
$\displaystyle+p(a=y,b=\bar{x}|x,y)P(\left|x\right\rangle\otimes\left|+\right\rangle\otimes\left|-\right\rangle\otimes\left|y\right\rangle)$
$\displaystyle+p(a=\bar{y},b=x|x,y)P(\left|x\right\rangle\otimes\left|-\right\rangle\otimes\left|+\right\rangle\otimes\left|y\right\rangle)$
$\displaystyle+p(a=\bar{y},b=\bar{x}|x,y)P(\left|x\right\rangle\otimes\left|+\right\rangle\otimes\left|+\right\rangle\otimes\left|y\right\rangle)\big{)}\,,$

where for brevity, given a state $\left|\psi\right\rangle$ we have defined the
projector
$P(\left|\psi\right\rangle)=\left|\psi\right\rangle\left\langle\psi\right|$.
Thus, by noting that $\text{Tr}\left\\{H_{q}\left|x\right\rangle\left\langle
x\right|\right\\}=\epsilon/2$, with $x=0,1$ and
$\text{Tr}\left\\{H_{q}\left|\pm\right\rangle\left\langle\pm\right|\right\\}=(1\pm
1)\epsilon/2$, from $\langle E\rangle=\text{Tr}\left\\{\rho H\right\\}$ we get
Eq. (13). Similarly, we deduce that the final energy of the circle qubits has
the probability distribution

$p(E_{c})=\sum_{i=0}^{2}p_{i}\delta(E_{c}-i\epsilon)\,,$ (15)

where $p_{0}=p_{succ}$. For a definite causal order, the optimal process
described above gives $p_{succ}=1/2$, and no chance of getting both bits
wrong, thus $p_{2}=0$ and $p_{1}=1/2$. In this case, the energy of the circle
qubits is $\epsilon$ or zero with the same probability. For an indefinite
causal order such that $p_{succ}>1/2$, the probability that this energy is
zero is larger than any definite causal order. Our next question is if it is
possible to achieve a gain in the average extracted work, i.e., to get
$\langle E\rangle<3\epsilon/2$. Surprisingly, for the system under
consideration, the answer is negative, because of the presence of a non-zero
probability $p_{2}$. In particular, we find the upper bound (see Appendix A
for the proof)

$p_{succ}-p_{2}\leq 1/2\,,$ (16)

from which $\langle E\rangle\geq 3\epsilon/2$. Thus, in the presence of
indefinite causal order, if $p_{succ}>1/2$ then $p_{2}>0$ and $p_{1}<1/2$.
Concerning the energy of the circle qubits, with respect to the optimal
causally separable process, we get a larger probability $p_{succ}$ that is
zero, a smaller probability $p_{1}$ that is $\epsilon$, but the non-zero
probability $p_{2}$ that is $2\epsilon$ gives an average extracted work
$\langle w\rangle=2\epsilon-\langle E\rangle$ not larger than $\epsilon/2$,
and so also a larger variance. We note that this result is quite general,
i.e., it is not possible to get a gain in the average work extracted. In
general, by performing a measurement of the spin with respect to an arbitrary
direction, the average energy of the square qubits does not change, and the
energy of the circle qubits still has the probability distribution of Eq. (17)
for optimal local cyclic unitary transformations. However, if interactions
between the parties are allowed, in particular if we consider the interaction
between the circle qubits
$H_{int}=-\epsilon\left|+\right\rangle\left\langle+\right|\otimes\left|+\right\rangle\left\langle+\right|$,
then the final energy of the circle qubits has the probability distribution

$p(E_{c})=p_{succ}\delta(E_{c})+(1-p_{succ})\delta(E_{c}-\epsilon)\,,$ (17)

which depends only on the probability $p_{succ}$. In this case, since the
average energy of the square qubits does not change, by considering that the
interaction lowers the initial average energy to $2\epsilon-\epsilon/4$, the
average work extracted is $\langle w\rangle=(p_{succ}-1/4)\epsilon$, which of
course, in the presence of an indefinite causal order such that
$p_{succ}>1/2$, is larger than the causal bound $\epsilon/4$. We note that if
the qubits are realized by using spinless fermions, the Hamiltonian reads
$H=H_{red}+H_{blue}+H_{int}$, where $H_{\alpha}=\epsilon
n_{\alpha}^{A}+\epsilon n_{\alpha}^{B}$ and $H_{int}=-\epsilon
n_{blue}^{A}n_{red}^{B}$, where $n_{\alpha}^{A}$ and $n_{\alpha}^{B}$ are the
number operators of the fermions of Alice and Bob, respectively, with
$\alpha=red,blue$. Of course, to implement the local cyclic unitary
transformations the interaction $H_{int}$ is switched off in the corresponding
time interval, and switched on at the end.

In the end, it is worth observing that the work extraction can be related to
the formation of correlations between the red qubits and the blue ones in the
final state $\rho$ of Eq. (12). In particular, for the optimal process with
definite causal order giving $p_{succ}=1/2$, in the final state the red qubits
are not correlated with the blue ones. In general, the total correlations are
quantified by the mutual information
$I_{red:blue}=S_{red}+S_{blue}-S_{red,blue}$, where $S_{red}$ and $S_{blue}$
are the von-Neumann entropies of the reduced final states of the red and blue
qubits, respectively. In detail, the von-Neumann entropy of a state $\rho$ is
defined as $S=-\text{Tr}\left\\{\rho\log_{2}\rho\right\\}$. Conversely,
$S_{red,blue}$ is the von-Neumann entropy of the final state $\rho$ of the
total system, which explicitly reads
$S_{red,blue}=-\sum_{a,b,x,y}\lambda_{abxy}\log_{2}\lambda_{abxy}$, where
$\lambda_{abxy}=p(x,y)p(a,b|x,y)$ are the eigenvalues of the density matrix
$\rho$. Since $\sum_{a,b}\lambda_{abxy}=p(x,y)=1/4$ and
$\sum_{x,y}\lambda_{abxy}=p^{\prime}(a,b)$, where $p^{\prime}(a,b)$ is the
probability to get the outputs $a$ and $b$, we get
$S_{red,blue}=2+H_{A,B}-I_{I:O}$, where
$H_{A,B}=-\sum_{a,b}p^{\prime}(a,b)\log_{2}p^{\prime}(a,b)$ is the Shannon
entropy corresponding to the outputs $a$, $b$, and $I_{I:O}$ is the mutual
information between the inputs $x$,$y$ and the outputs $a$, $b$. Then, we get

$I_{red:blue}=S_{red}+S_{blue}-2+I_{I:O}-H_{A,B}\,,$ (18)

from which, since $S_{red}\leq 2$ and $S_{blue}\leq 2$, we get the upper bound

$I_{red:blue}\leq 2+I_{I:O}-H_{A,B}\,,$ (19)

where $0\leq H_{A,B}-I_{I:O}\leq 2$. We note that Eq. (19) connects the
formation of correlations between the red and blue qubits only to the
information exchange. Concerning $H_{A,B}$, as shown in Ref. Francica22 ,
there exist causally non-separable process matrices $W$ such that
$H_{A,B}(W)>\max_{W_{sep}\in\mathcal{S}_{W}}H_{A,B}(W_{sep})$, where
$\mathcal{S}_{W}$ is the set of all the separable processes $W_{sep}$, defined
by Eq. (11), such that $\Delta(W^{A\prec B})=\Delta(W)$ and $\Delta(W^{B\prec
A})=\Delta(W)$, and $\Delta(W)$ is the non-signalling part of the process
matrix $W$ defined as $\Delta(W)={{}_{A_{O}B_{O}}W}$. However, this does not
rule out that $H_{A,B}-I_{I:O}$ is smaller for an indefinite causal order,
getting a larger amount of correlations between the red and the blue qubits.
For instance, we consider the process matrix

$W=\frac{1}{4}\left(I^{\otimes
4}+\alpha(\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\sigma_{z}^{B_{I}}I^{B_{O}}+\sigma_{z}^{A_{I}}I^{A_{O}}\sigma_{x}^{B_{I}}\sigma_{x}^{B_{O}})\right)\,,$
(20)

with $0\leq\alpha\leq 1/\sqrt{2}$, where the tensor products are implicit, and
the local operations

$\displaystyle M^{A_{I}A_{O}}_{0|0}$ $\displaystyle=$ $\displaystyle
M^{B_{I}B_{O}}_{0|0}=0\,,$ (21) $\displaystyle M^{A_{I}A_{O}}_{1|0}$
$\displaystyle=$ $\displaystyle
M^{B_{I}B_{O}}_{1|0}=\left|\varphi^{+}\right\rangle\left\langle\varphi^{+}\right|\,,$
(22) $\displaystyle M^{A_{I}A_{O}}_{0|1}$ $\displaystyle=$ $\displaystyle
M^{B_{I}B_{O}}_{0|1}=\left|0\right\rangle\left\langle
0\right|\otimes\left|0\right\rangle\left\langle 0\right|\,,$ (23)
$\displaystyle M^{A_{I}A_{O}}_{1|1}$ $\displaystyle=$ $\displaystyle
M^{B_{I}B_{O}}_{1|1}=\left|1\right\rangle\left\langle
1\right|\otimes\left|0\right\rangle\left\langle 0\right|\,.$ (24)

We get $p_{succ}=5(1+\alpha)/16$, thus the violation of causal inequalities
becomes larger as $\alpha$ increases, for $\alpha=1/\sqrt{2}$ is maximum and
we get the causally non-separable process matrix of Ref. Branciard16 . We find
that $H_{A,B}$ increases, $H_{A,B}-I_{I:O}$ decreases and $I_{red:blue}$
increases as $\alpha$ increases. In particular, for $\alpha=1/\sqrt{2}$, we
get a non-zero mutual information $I_{red:blue}\approx 1.0951$. Moreover, the
bound is $I_{I:O}-H_{A,B}+2\approx 1.2993$, so that $I_{red:blue}$ is close to
the upper bound since $S_{red}=S_{blue}\approx 1.8979$ is close to two.
Furthermore, we get $p_{2}=(5+\alpha)/16$, and so $p_{succ}-p_{2}=\alpha/4$,
in agreement with Eq. (16). We note that the bound of Eq. (16) can be
saturated by choosing opportunely the local operations, e.g., for
$\alpha=1/\sqrt{2}$ we get $p_{succ}=1/2$ and $p_{2}=0$ for the operations

$\displaystyle M^{A_{I}A_{O}}_{0|0}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}+\sigma_{z}^{A_{I}}I^{A_{O}}+I^{A_{I}}\sigma_{z}^{A_{O}}+\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\right)\,,$
(25) $\displaystyle M^{A_{I}A_{O}}_{1|0}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}-\sigma_{z}^{A_{I}}I^{A_{O}}-I^{A_{I}}\sigma_{z}^{A_{O}}+\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\right)\,,$
(26) $\displaystyle M^{A_{I}A_{O}}_{0|1}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}+\sigma_{z}^{A_{I}}I^{A_{O}}-I^{A_{I}}\sigma_{z}^{A_{O}}-\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\right)\,,$
(27) $\displaystyle M^{A_{I}A_{O}}_{1|1}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}-\sigma_{z}^{A_{I}}I^{A_{O}}+I^{A_{I}}\sigma_{z}^{A_{O}}-\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\right)\,,$
(28) $\displaystyle M^{B_{I}B_{O}}_{0|0}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}+\frac{1}{\sqrt{2}}M^{B_{I}B_{O}}_{+++}+I^{B_{I}}\sigma_{x}^{B_{O}}\right)\,,$
(29) $\displaystyle M^{B_{I}B_{O}}_{1|0}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}-\frac{1}{\sqrt{2}}M^{B_{I}B_{O}}_{-+-}+I^{B_{I}}\sigma_{x}^{B_{O}}\right)\,,$
(30) $\displaystyle M^{B_{I}B_{O}}_{0|1}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}+\frac{1}{\sqrt{2}}M^{B_{I}B_{O}}_{+--}-I^{B_{I}}\sigma_{x}^{B_{O}}\right)\,,$
(31) $\displaystyle M^{B_{I}B_{O}}_{1|1}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(I^{\otimes
2}-\frac{1}{\sqrt{2}}M^{B_{I}B_{O}}_{+++}+I^{B_{I}}\sigma_{x}^{B_{O}}\right)\,,$
(32)

where we have defined
$M^{B_{I}B_{O}}_{\pm\pm\pm}=\sigma_{z}^{B_{I}}I^{B_{O}}\pm\sigma_{x}^{B_{I}}I^{B_{O}}\pm\sigma_{z}^{B_{I}}\sigma_{x}^{B_{O}}\pm\sigma_{x}^{B_{I}}\sigma_{x}^{B_{O}}$.

##  IV Conclusions

To summarize, we provided a general scheme with the aim to achieve an
advantage in work extraction games in the presence of indefinite causal order.
For a system of four qubits, we showed that the probability that the circle
qubits have a zero final energy in a single realization is equal to the
success probability of the “guess your neighbor’s input” game. Anyway, for our
scheme, the average work extracted can not be larger than the one extracted
for definite causal order if the qubits are not interacting. This suggests
that, in this case, non-causal process matrices do not give a violation of the
second law of thermodynamics corresponding to a definite causal order, which
in general bounds the amount of extractable work. However, it is possible to
add an interaction between the qubits to get a larger average work with
respect to any causal process matrix. Anyway, we note that if interactions are
allowed, for a definite causal order in principle Bob (or Alice) can realize a
global operation to always lower to zero the energy of the circle qubits.
Concerning the final state obtained, the scheme leads to the formation of
correlations between the two couples of the two qubits in the presence of
indefinite causal order. In conclusion, we hope that our results can open a
new avenue in applying indefinite causal order structures and can inspire
further investigations and applications. In particular, further investigations
are needed to understand how the thermodynamic time arrow given by the second
law can be affected by the causal structure.

### Acknowledgements

The author acknowledges financial support from the project BIRD 2021
”Correlations, dynamics and topology in long-range quantum systems” of the
Department of Physics and Astronomy, University of Padova.

##  Appendix A Upper bound

To prove the bound of Eq. (16), we note that by using Eq. (1) we get

$\displaystyle p_{succ}-p_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\text{Tr}\bigg{\\{}\bigg{(}\left(M^{A_{I}A_{O}}_{0|0}-M^{A_{I}A_{O}}_{1|1}\right)\left(M^{B_{I}B_{O}}_{0|0}-M^{B_{I}B_{O}}_{1|1}\right)$
(33)
$\displaystyle-\left(M^{A_{I}A_{O}}_{0|1}-M^{A_{I}A_{O}}_{1|0}\right)\left(M^{B_{I}B_{O}}_{0|1}-M^{B_{I}B_{O}}_{1|0}\right)\bigg{)}W\bigg{\\}}\,,$

where the tensor products are implicit, from which

$\displaystyle p_{succ}-p_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\text{Tr}\bigg{\\{}\bigg{(}\left(M^{A_{I}A_{O}}_{0}-M^{A_{I}A_{O}}_{1}\right)\left(M^{B_{I}B_{O}}_{0|1}-M^{B_{I}B_{O}}_{1|0}\right)$
(34)
$\displaystyle+\left(M^{A_{I}A_{O}}_{0|1}-M^{A_{I}A_{O}}_{1|0}\right)\left(M^{B_{I}B_{O}}_{0}-M^{B_{I}B_{O}}_{1}\right)\bigg{)}W\bigg{\\}}\,,$

where $M^{A_{I}A_{O}}_{x}=\sum_{a}M^{A_{I}A_{O}}_{a|x}$ and
$M^{B_{I}B_{O}}_{y}=\sum_{b}M^{B_{I}B_{O}}_{b|y}$, which are such that
$\text{Tr}\left\\{M^{A_{I}A_{O}}_{x}M^{B_{I}B_{O}}_{y}W\right\\}=1$ for any
$x$ and $y$ (see, e.g., Ref. Oreshkov12 ). In general, the process matrix $W$
can be written as

$W=\Delta(W)+\sum_{\alpha,i,j}c_{\alpha
ij}\sigma^{A_{I}}_{\alpha}\sigma^{A_{O}}_{i}\sigma^{B_{I}}_{j}I^{B_{O}}+c^{\prime}_{i\alpha
j}\sigma^{A_{I}}_{i}I^{A_{O}}\sigma^{B_{I}}_{\alpha}\sigma^{B_{O}}_{j}\,,$
(35)

where $\Delta(W)={{}_{A_{O}B_{O}}W}$ and $c_{\alpha ij}$ and
$c^{\prime}_{i\alpha j}$ are real parameters. We have that
$\text{Tr}_{A_{O}}\left\\{O^{A_{I}}I^{A_{O}}M^{A_{I}A_{O}}_{x}\right\\}=O^{A_{I}}$
for any $O^{A_{I}}$, and a similar equation for $M^{B_{I}B_{O}}_{y}$, so that

$p_{succ}-p_{2}=\frac{1}{4}\sum_{\alpha,i,j}c_{\alpha ij}m_{\alpha
ij}+c^{\prime}_{i\alpha j}m^{\prime}_{i\alpha j}\,,$ (36)

where we have defined

$\displaystyle m_{\alpha ij}$ $\displaystyle=$
$\displaystyle\text{Tr}\left\\{\left(M^{A_{I}A_{O}}_{0}-M^{A_{I}A_{O}}_{1}\right)\sigma^{A_{I}}_{\alpha}\sigma^{A_{O}}_{i}\right\\}$
(37)
$\displaystyle\times\text{Tr}\left\\{\left(M^{B_{I}B_{O}}_{0|1}-M^{B_{I}B_{O}}_{1|0}\right)\sigma^{B_{I}}_{j}I^{B_{O}}\right\\}\,,$
$\displaystyle m^{\prime}_{i\alpha j}$ $\displaystyle=$
$\displaystyle\text{Tr}\left\\{\left(M^{A_{I}A_{O}}_{0|1}-M^{A_{I}A_{O}}_{1|0}\right)\sigma^{A_{I}}_{i}I^{A_{O}}\right\\}$
(38)
$\displaystyle\times\text{Tr}\left\\{\left(M^{B_{I}B_{O}}_{0}-M^{B_{I}B_{O}}_{1}\right)\sigma^{B_{I}}_{\alpha}\sigma^{B_{O}}_{j}\right\\}\,.$

We note that the operator $M^{A_{I}A_{O}}_{a|x}$ can be expressed as

$M^{A_{I}A_{O}}_{a|x}=\frac{q_{a|x}}{2}\bigg{(}I^{\otimes
2}+\vec{r}_{a|x}\cdot\vec{\sigma}^{A_{I}}I^{A_{O}}+I^{A_{I}}\vec{\sigma}^{A_{O}}\cdot\vec{s}_{a|x}+\sum_{i,j}t^{a|x}_{ij}\sigma_{i}^{A_{I}}\sigma_{j}^{A_{O}}\bigg{)}\,,$
(39)

where $0\leq q_{a|x}\leq 1$, $\sum_{a}q_{a|x}=1$ and
$\sum_{a}q_{a|x}\vec{r}_{a|x}=0$ since
$\text{Tr}_{A_{O}}\left\\{M^{A_{I}A_{O}}_{x}\right\\}=I^{A_{I}}$. Thus, we get

$\displaystyle\sum_{\alpha,i}\left|\text{Tr}\left\\{\left(M^{A_{I}A_{O}}_{0}-M^{A_{I}A_{O}}_{1}\right)\sigma^{A_{I}}_{\alpha}\sigma^{A_{O}}_{i}\right\\}\right|^{2}$
$\displaystyle=$ $\displaystyle
4\left(||\vec{s}_{0}-\vec{s}_{1}||^{2}+\sum_{i,j}\left|t^{0}_{ij}-t^{1}_{ij}\right|^{2}\right)$
(40) $\displaystyle\leq$ $\displaystyle 16\,,$ (41)

where $\vec{s}_{x}=\sum_{a}q_{a|x}\vec{s}_{a|x}$ and
$t^{x}_{ij}=\sum_{a}q_{a|x}t^{a|x}_{ij}$, and

$\displaystyle\sum_{i}\left|\text{Tr}\left\\{\left(M^{A_{I}A_{O}}_{0|1}-M^{A_{I}A_{O}}_{1|0}\right)\sigma^{A_{I}}_{i}I^{A_{O}}\right\\}\right|^{2}$
$\displaystyle=$ $\displaystyle
4||q_{0|1}\vec{r}_{0|1}-q_{1|0}\vec{r}_{1|0}||^{2}$ (42) $\displaystyle\leq$
$\displaystyle 4\,.$ (43)

Thus, we define the vectors $\vec{m}$ and $\vec{m}^{\prime}$ with components
$m_{\alpha ij}$ and $m^{\prime}_{i\alpha j}$, respectively, and we get
$||\vec{m}||\leq 8$ and $||\vec{m}^{\prime}||\leq 8$. From Eq. (36), we get

$|p_{succ}-p_{2}|=\frac{1}{4}\left|\vec{c}\cdot\vec{m}+\vec{c}^{\prime}\cdot\vec{m}^{\prime}\right|\leq\frac{1}{4}\left(|\vec{c}\cdot\vec{m}|+|\vec{c}^{\prime}\cdot\vec{m}^{\prime}|\right)\,,$
(44)

where we have defined the vectors $\vec{c}$ and $\vec{c}^{\prime}$ with
components $c_{\alpha ij}$ and $c^{\prime}_{i\alpha j}$, respectively. To find
an upper bound of
$|\vec{c}\cdot\vec{m}|+|\vec{c}^{\prime}\cdot\vec{m}^{\prime}|$, we consider
the case where $\vec{c}$ has only one non-zero component, which is $c/4$, and
$\vec{c}^{\prime}$ has only one non-zero component, which is $c^{\prime}/4$.
Since $W\geq 0$, there are two cases: or $|c|\leq|\cos\theta|$ and
$|c^{\prime}|\leq|\sin\theta|$ for a certain $\theta$, e.g., for process
matrices as

$W=\frac{1}{4}\left(I^{\otimes
4}+c\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\sigma_{z}^{B_{I}}I^{B_{O}}+c^{\prime}\sigma_{z}^{A_{I}}I^{A_{O}}\sigma_{x}^{B_{I}}\sigma_{x}^{B_{O}}\right)\,,$
(45)

or $|c|+|c^{\prime}|\leq 1$, e.g., for process matrices as

$W=\frac{1}{4}\left(I^{\otimes
4}+c\sigma_{z}^{A_{I}}\sigma_{z}^{A_{O}}\sigma_{z}^{B_{I}}I^{B_{O}}+c^{\prime}\sigma_{z}^{A_{I}}I^{A_{O}}\sigma_{z}^{B_{I}}\sigma_{z}^{B_{O}}\right)\,.$
(46)

It is easy to see that in both cases,
$|\vec{c}\cdot\vec{m}|+|\vec{c}^{\prime}\cdot\vec{m}^{\prime}|\leq 2$, from
which $|p_{succ}-p_{2}|\leq 1/2$.

## References

  * (1) O. Oreshkov, F. Costa and Č. Brukner, Nat. Commun. 3, 1092 (2012).
  * (2) C. Branciard, M. Araújo, A. Feix, F. Costa and Č. Brukner, New J. Phys. 18, 013008 (2016).
  * (3) A. Feix, M. Araújo and Č. Brukner, New J. Phys. 18, 083040 (2016).
  * (4) G. Chiribella, G. M. D’Ariano, P. Perinotti and B. Valiron, Phys. Rev. A 88, 022318 (2013).
  * (5) E. Adlam, arXiv:2208.02721.
  * (6) A. Feix, M. Araújo and Č. Brukner, Phys. Rev. A 92, 052326 (2015).
  * (7) P. A. Guérin, A. Feix, M. Araújo and Č. Brukner, Phys. Rev. Lett. 117, 100502 (2016).
  * (8) D. Ebler, S. Salek and G. Chiribella, Phys. Rev. Lett. 120, 120502 (2018).
  * (9) G. Chiribella, Phys. Rev. A 86, 040301(R) (2012).
  * (10) M, Araújo, F. Costa and Č. Brukner, Phys. Rev. Lett. 113, 250402 (2014).
  * (11) L. M. Procopio, A. Moqanaki, M. Araújo, F. Costa, I. Alonso Calafell, E. G. Dowd, D. R. Hamel, L. A. Rozema, Č. Brukner and P. Walther, Nat. Commun. 6, 7913 (2015).
  * (12) G. Rubino, L. A. Rozema, A. Feix, M. Araújo, J. M. Zeuner, L. M. Procopio, Č. Brukner and P. Walther, Sci. Adv. 3, e1602589 (2017).
  * (13) K. Goswami, C. Giarmatzi, M. Kewming, F. Costa, C. Branciard, J. Romero and A. G. White, Phys. Rev. Lett. 121, 090503 (2018).
  * (14) K. Goswami, Y. Cao, G. A. Paz-Silva, J. Romero and A. G. White, Phys. Rev. Research 2, 033292 (2020).
  * (15) G. Rubino, L. A. Rozema, D. Ebler, H. Kristjánsson, S. Salek, P. A. Guérin, A. A. Abbo, C. Branciard, Č. Brukner, G. Chiribella and P. Walther, Phys. Rev. Research 3, 013093 (2021).
  * (16) X. Nie, X. Zhu, C. Xi, X. Long, Z. Lin, Y. Tian, C. Qiu, X. Yang, Y. Dong, J. Li, T. Xin, and D. Lu, arXiv:2011.12580.
  * (17) H. Cao, N. Wang, Z.-A. Jia, C. Zhang, Y. Guo, B.-H. Liu, Y. F. Huang, C.-F. Li, and G.-C. Guo, arXiv:2101.07979.
  * (18) D. Felce, V. Vedral, and F. Tennie, arXiv:2107.12413.
  * (19) D. Felce, and V. Vedral, Phys. Rev. Lett. 125, 070603 (2020).
  * (20) T. Guha, M. Alimuddin, and P. Parashar, Phys. Rev. A 102, 032215 (2020).
  * (21) K. Simonov, G. Francica, G. Guarnieri, and M. Paternostro, Phys. Rev. A 105, 032217 (2022).
  * (22) P. R. Dieguez, V. F. Lisboa, and R. M. Serra, arXiv:2205.14406
  * (23) M. Capela, H. Verma, F. Costa, and L. C. Céleri, arXiv:2208.03205.
  * (24) T. Guha, S. Roy, K. Simonov, and Z. Zimborás, arXiv:2208.04034.
  * (25) G. Francica, J. Goold, F. Plastina, and M. Paternostro, npj Quantum Inf. 3, 12 (2017).
  * (26) A. E. Allahverdyan, R. Balian and Th.M. Nieuwenhuizen, Europhys. Lett. 67, 4 (2004).
  * (27) M.D. Choi, Linear Algebra Appl. 10, 285-90 (1975).
  * (28) A. Jamiołkowski, Rep. Math. Phys. 3, 275-8 (1972).
  * (29) M. Araújo, C. Branciard, F. Costa, A. Feix, C. Giarmatzi and Č. Brukner, New J. Phys. 17, 102001 (2015).
  * (30) A. Baumeler, A. Feix and S. Wolf, Phys. Rev. A 90, 042106 (2014).
  * (31) O. Oreshkov and C. Giarmatzi, New J. Phys. 18, 093020 (2016).
  * (32) G. Francica, Phys. Lett. A 422 (2022) 127815; G. Francica, Phys. Lett. A 436 (2022) 128087.


#  DIAS-STP-24-01 UCI-HEP-TR-2024-03
Constraints on Variation of the Weak Scale
from Big Bang Nucleosynthesis

Anne-Katherine Burns<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of California, Irvine, CA 92697 USA  Venus Keus<EMAIL_ADDRESS>Dublin Institute for Advanced Studies, School of Theoretical Physics, 10
Burlington road, Dublin, D04 C932, Ireland  Marc Sher<EMAIL_ADDRESS>High
Energy Theory Group, Physics Department, William & Mary, Williamsburg, VA
23187, USA  Tim M.P. Tait<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of California, Irvine, CA 92697 USA

###### Abstract

Recently, the EMPRESS collaboration has included new data in the extraction of
the primordial 4He abundance from Big Bang Nucleosynthesis (BBN), resulting in
a determination that differs from the previous value and from theoretical
expectations. There have been several studies attempting to explain this
anomaly which involve variation of fundamental constants between the time of
BBN and the present. Since the Higgs vacuum expectation value (vev) is the
only dimensionful parameter of the Standard Model and it is already known to
vary during the electroweak phase transition, we consider the possibility that
the vev is slightly different during BBN compared to its present value. A
modification of the vev changes not only particle masses but also affects,
through mass thresholds, the QCD confinement scale. We use the recently
developed PRyMordial program to study this variation and its impact on the 4He
and deuterium abundances. We find that bounds on $|\frac{\delta v}{v}|$ are
approximately $0.01$, and that the EMPRESS result can be explained within
$2\sigma$ if $-0.03<\frac{\delta v}{v}<-0.01$, but at the cost of worsening
the current $2\sigma$ discrepancy in the deuterium abundance to over
$3\sigma$.

##  I Introduction

The $\Lambda CDM$ model has been extraordinarily successful in describing the
cosmological history and evolution of the Universe. However there are some
anomalies in the model. The most well-known is the Hubble tension, in which
the value of the Hubble constant measured from the cosmic microwave background
(CMB) by Planck Planck:2018vyg , $H_{0}=67.4\pm 0.5\ {\rm km/s/Mpc}$ differs
by $5\sigma$ from the measurement of Cepheids and type Ia supernovae by the
SH0ES collaborationRiess:2021jrx , $H_{0}=73.0\pm 1.0\ {\rm km/s/Mpc}$. There
have been numerous studies trying to explain this discrepancy (see Ref.
Dainotti:2023yrk  for an extensive review). Another anomaly concerns the 4He
abundance, $Y_{p}$ from Big Bang Nucleosynthesis (BBN). A recent report by
EMPRESS Matsumoto:2022tlr  adds measurements from 10 additional extremely
metal-poor galaxies which, combined with the previously existing dataset,
results in a fit to $Y_{p}=0.2370\pm 0.0033$ that differs from the value
obtained based only the pre-existing data Kurichin:2021ppm ; Hsyu:2020uqb .
This has led to several studies, such as the implications of the result for
the lepton asymmetry Kawasaki:2022hvx ; Burns:2022hkq ; Borah:2022uos ;
Escudero:2022okz

The variation of fundamental constants is a subject that has been around since
Dirac’s large number hypothesis dirac . Most work in the past has focused on
variation in the fine-structure constant Bergstrom:1999wm ; Nollett:2002da ;
Mosquera:2007vr  and the electron mass Barrow:2005qf ; Scoccola:2008yf  and
both have been used to study the above anomalies. Variation of the fine-
structure constant to solve the Hubble tension was proposed in a mirror dark
sector model recently Zhang:2022ujw . Variation of the electron mass
Barrow:2005qf  has also been proposed as a solution to the Hubble tension in
Ref. Seto:2022xgx  in which they also studied the effects on Big Bang
Nucleosynthesis; a mechanism for such a variation was proposed in Ref.
Solomon:2022qqf . A detailed series of papers studying variation of both the
fine-structure constant and the electron mass is by Hart and Chluba
Hart:2017ndk ; Hart:2019dxi ; Hart:2021kad . The helium anomaly is more
recent. It was shown by Seto, Takahashi and Toda Seto:2023yal  that a
variation in the fine-structure constant alone would be sufficient to explain
the anomaly.

While one can consider variation of the fine-structure constant and the
electron mass, we feel that it would be more reasonable to consider variation
in the vacuum expectation value (vev) of the Higgs field (or, equivalently,
the Higgs mass term). The vev is the only dimensionful parameter in the
Standard Model, and it is also the most mysterious, in the sense that its
value is many, many orders of magnitude smaller that one would naturally
expect. In addition, we already know that the Standard Model vev is
cosmologically dynamical, and varies substantially during the electroweak
phase transition. Thus it would seem to be more reasonable to consider
variations in the vev rather than in one of the many dimensionless parameters.
One can imagine an extended model in which the Higgs vev slowly rolls to its
current value but has not quite reached it at the MeV temperature scale.

The possibility effects of a changing vev on BBN was first discussed in 1988
by Dixit and Sher Dixit:1987at , but this was a crude calculation that only
considered the helium abundance. There have since been several studies
Yoo:2002vw ; Chamoun:2005xr ; Gassner:2006fr ; Landau:2008re  of BBN due to a
changing vev and others Muller:2004gu ; Coc:2006sx  on the effects of changing
many parameters, including the vev. These papers all include the effects not
only on the electron mass but also on quark masses, which then affects the
pion and other meson masses and thus the strong force (leading, for example,
to changes in the deuteron binding energy). Other works did not directly
discuss modifying the vev, but only $\Lambda_{QCD}$ Kneller:2003xf , the
deuteron binding energy Dmitriev:2003qq  and quark masses Bedaque:2010hr .
Like these works, our approach is phenomenological and we will not present a
specific model with a changing vev at the time of BBN. A recent paper
Laverda:2024qjt  did present a model in which the Higgs field may not achieve
a thermal spectrum at the time of nucleosynthesis - the effects on the vev are
unclear. Further investigation of this model would be interesting.

However, these works did not include the effects of a changing Higgs vev on
the value of $\Lambda_{QCD}$, which clearly will affect strong interaction
dynamics. As noted by Agrawal, et al. Agrawal:1997gf ; Agrawal:1998xa , a
variation in the Higgs vev will modify all of the quark masses and this in
turn will impact the running of the strong coupling constant through the quark
thresholds. The effect is somewhat smaller (they find $\Lambda_{QCD}$ varying
by roughly $(v/v_{0})^{0.25})$), but can still be substantial. These effects
were included in Ref. Jeltema:1999na  which was an anthropic principle based
study of the triple alpha process in stars, but have not been included in
analyses of BBN.

Since these above works, we have learned much more about the current baryon-
photon ratio, the nuclear abundances and have better codes Burns:2023sgx . One
such code, PRyMordial was created to fill several gaps that currently exist in
the landscape of publicly available BBN codes. Namely, the code allows for
easy and flexible exploration of a variety of new physics models, including
the investigation of varying input parameters on both $N_{eff}$ and the final
abundance values. Prymordial also allows users to examine the uncertainties in
different sets of thermonuclear reaction rates and to scale these chosen
reaction rates as desired. In addition, PRyMordial calculates the
thermodynamics of the plasma from first principles, the results of which serve
as the initial conditions for the full calculation of neutron to proton
conversion. Unlike other codes, PRyMordial is written in python, allowing its
users to perform Monte Carlo (MC) analyses with ease by interfacing with
standard MC libraries. In order to maximize efficiency, the code also has an
option to use Julia to speed up the computation. PRyMordial has already been
used by several groups to study the effects of a variety of new physics
scenarios on BBN Burns:2022hkq ; Meissner:2023voo ; Schoneberg:2024ifp ;
Giri:2023pyy ; Hong:2023fcy ; Bianchini:2023ubu ; Chowdhury:2022ahn . In this
paper, using PRyMordial, we will calculate the abundances of light nuclei as a
function of the Higgs vev including all of these effects and in light of the
EMPRESS result.

##  II Procedure

When the vev is varied, the quark and electron masses receive fractional
corrections of ${\delta v}/{v}$. While the bulk of the neutron and proton
masses are the result of confinement, their masses, and most importantly their
mass difference, is sensitive to the change in quark masses. As is well known,
the down quark mass is roughly $2.5$ MeV heavier than the up quark, but the
neutron is only $1.293$ MeV heavier than the proton. The discrepancy is due to
the electromagnetic energy of the up quarks within the proton, which is thus
approximately $1.2$ MeV, independent of the Higgs vev. As a result, we take
the neutron-proton mass difference in MeV to be

$\frac{m_{n}-m_{p}}{\rm MeV}=2.493~{}(1+\frac{\delta v}{v})-1.2$ (1)

Weinberg Weinberg:1972kfs  lists the various reaction rates for the six
processes involving $n,p,e,\nu$ in nucleosynthesis and one can see that the
variation in the neutron-proton mass difference and the electron mass will
suffice to determine the variation in these processes, including, of course,
the neutron lifetime. However this was a fairly crude calculation of the 4He
abundance - to get a more precise calculation and to include the deuteron,
triton and 3He, one needs to look at nuclear processes. The binding energy of
the deuteron is especially significant since this leads to the well-known
deuterium bottleneck.

The pion mass is particularly important for the deuteron binding energy. The
pion mass in QCD is given by $m_{\pi}^{2}\simeq(m_{u}+m_{d})f_{\pi}$. While
the quark masses scalar linearly with $\delta v$, the value of $f_{\pi}$ is
proportional to $\Lambda_{QCD}$. As noted earlier, this will scale differently
from the quark masses and will depend on the vev through mass thresholds. In
Appendix A, we show that $\Lambda_{QCD}$ scales as $(\delta v)^{0.25}$ and
thus the pion mass scales as $(\delta v)^{1.25/2}$.

We now turn to the deuteron binding energy. Since we are interested in the
impact of a small change in the vev, high precision in the standard
calculation is not necessary. We model the nucleon-nucleon potential for the
deuteron as an exchange of pion, $\sigma$ and $\omega$ mesons as discussed in
the review article by Meissner Meissner:1987ge . The pion mass, as discussed
above, scales as $\left({\delta v}/{v}\right)^{1.25/2}$. The $\omega$ mass
comes primarily from QCD and so scales as $\left({\delta
v}/{v}\right)^{0.25}$. The $\sigma$ is believed to be a two-pion correlated
state. Lin and Serot Lin:1990cx  calculated the $\sigma$ mass in terms of the
pion-nucleon coupling, the pion mass and the nucleon mass. Varying the masses
in their expressions, we find that the $\sigma$ mass is insensitive to the
pion mass, scaling as $\Lambda_{QCD}$. The combined exchange results in the
sum of three Yukawa potentials - repulsive from the $\omega$ and attractive
from $\sigma$ and the pion.

Of course, there are other nuclei which participate in nucleosynthesis. We
assume that the binding energies and matrix elements involving these heavier
nuclei scale as $\Lambda_{QCD}$, which is expected to be a decent
approximation since the pion interactions are long-range.

##  III Results

We begin with the six processes involving $n,p,e,\nu$. The expressions for the
total $n\to p$ and $p\to n$ conversion rates are given in Weinberg
Weinberg:1972kfs . Using Eq. (1) for the neutron-proton mass difference and
scaling the electron mass fractionally as ${\delta v}/{v}$, the relationship
between the electron and neutrino temperatures is slightly modified. Putting
these effects together with the factor of $1/(1+d)^{4}$ coming from the shift
in the $W$ mass, the resulting fractional change for the $n\to p$ and $p\to n$
conversion rates as a function of $T_{\gamma}$, for $d\equiv|{\delta
v}/{v}|=0.01$, are shown in Figure 1. The fractional changes in the rates
scale approximately linearly with $\delta v$.

Figure 1: The fractional change in the $\lambda(n\to p)$ and $\lambda(p\to n)$
rates for $d\equiv|{\delta v}/{v}|=0.01$. Note that the fractional change in
$\lambda(p\to n)$ gets very large at low temperatures since the rate itself
vanishes at low temperature. For each rate, we include a factor of
$1/(1+d)^{4}$ coming from the shift in the $W$ mass.

As noted above, the deuteron binding energy also receives important
modifications from a change in the vev. There are two attractive Yukawa
potentials from the $\sigma$ and $\pi$ and a repulsive Yukawa potential from
the $\omega$. Each potential has a corresponding coupling constant indicated
by $g_{\pi}$, $g_{\sigma}$ and $g_{\omega}$. Families of viable parameters can
be found for $d=0$ by adjusting these three parameters such that the solution
to the Schrödinger equation gives a binding energy of $2.2$ MeV. One such
solution111We have examined other values of $g_{\pi}$, $g_{\sigma}$ and
$g_{\omega}$ which give the correct $d=0$ binding energy of $2.2$ MeV, and
find that they make very similar predictions for small $d$ to the ones we have
chosen., for example, has $g_{\pi}=11.97$, $g_{\sigma}=8.46$,
$g_{\omega}=21.19$. These are dimensionless and should not change much due to
small changes in the vev. Results for small $d\neq 0$ are obtained by
adjusting the meson masses as described above and recalculating the binding
energy. The resulting binding energies are presented in Figure 2 for positive
and negative values of ${\delta v}/{v}$. We do not consider values of
$|{\delta v}/{v}|>0.1$ which lead to either very large binding energies or an
unbounded deuteron for positive and negative values, respectively.

Figure 2: Value of the deuterium binding energy with varied Higgs vev for
positive and negative values of ${\delta v}/{v}$.

With the change in the $n\leftrightarrow p$ processes and the deuteron binding
energy in hand, we turn to the other nuclear processes in big bang
nucleosynthesis. There are 12 key reactions which are needed to predict
deuterium and helium abundances and an additional 51 reactions to accurately
predict lithium. These reaction rates are typically determined empirically,
making it difficult to predict their dependence on the vev from first
principles. However, they do not involve the weak interactions and are not as
sensitive to the pion mass and thus on dimensional grounds, we assume their
rates vary as $\Lambda_{QCD}$, or $v^{0.25}$.

We use PRyMordial to determine the abundances of primordial helium and
deuterium as a function of ${\delta v}/{v}$. The result for $Y_{p}$, which
characterizes the abundance of 4He, is shown in Figure 3 for both positive and
negative values of ${\delta v}/{v}$.

Figure 3: Value of the abundance of 4He, $Y_{p}$, with varying Higgs vev for
positive and negative values of ${\delta v}/{v}$. The red dashed line is the
standard result without a change in the Higgs vev and the blue and green lines
are the results with the change. The coral and gold boxes give the
experimental values from the Particle Data Group and from the recent EMPRESS
experiment.

The blue line corresponds to the predicted 4He abundance from PRyMordial. The
uncertainty for $Y_{p}$ is negligible. The experimental value from Refs.
Kurichin:2021ppm ; Hsyu:2020uqb  (as listed in the Particle Data Group
Workman:2022ynf  summary) and the $2\sigma$ uncertainties are shown shaded
green, and the recent result from EMPRESS Matsumoto:2022tlr  shaded blue. For
positive ${\delta v}/{v}$, the $2\sigma$ upper bound is approximately $0.003$,
whereas for negative ${\delta v}/{v}$, the magnitude must be less than $0.01$
based on the PDG result and less than $0.03$ using the EMPRESS result. If the
EMPRESS result turns out to be correct, $\delta v=0$ would be mildly excluded,
and $\delta v$ between $-0.01$ and $-0.03$ would be able to explain the 4He
abundance.

Figure 4: Value of the abundance of deuterium with varying Higgs vev for
positive and negative values of ${\delta v}/{v}$. The red dashed line is the
standard result without a change in the Higgs vev and the blue and green lines
are the results with the change. The coral box gives the experimental values
from the Particle Data Group. The error in the theoretical curves are similar
to the error in the experimental values, as discussed in the text.

In Figure 4, we show the results for the deuterium abundance. Here, the
theoretical calculation of the abundance has an uncertainty of $0.05$, which
is similar to the error in the experimental value. Thus the $\delta v=0$ limit
has roughly a $2\sigma$ discrepancy between the theoretical and experimental
values. For positive values of ${\delta v}/{v}$, we see that ${\delta v}/{v}$
must be less than approximately $0.06$, but for negative ${\delta v}/{v}$ the
magnitude must be substantially smaller. This is somewhat unfortunate since
negative ${\delta v}/{v}$ is needed to explain the EMPRESS result for $Y_{p}$,
but that would make the discrepancy with deuterium much worse. Thus changing
the vev can explain one, but not both, of the anomalies.

##  IV Conclusions

Originally motivated by anthropic arguments, many have studied the possibility
that some of the constants of nature are time-dependent. Astrophysical studies
have bounded some parameters as have some cosmological studies. Recently, the
EMPRESS collaboration reported a measurement of the primordial 4He abundance
which is somewhat over $3\sigma$ from the theoretical value. There is also a
small $2\sigma$ discrepancy in the deuterium abundance. In this paper, we have
considered the possibility that a difference in the Higgs vev between the time
of BBN and the present could be responsible for one or both of these
discrepancies. We choose the Higgs vev to vary because it is the only
dimensionful parameter of the Standard Model and is already known to vary
during the earlier electroweak phase transition.

This was studied some time ago, but previous authors did not include the fact
that the QCD scale will also vary if the quark masses vary, due to threshold
effects. They also did not use the latest BBN codes, which are more accurate
than previous ones. The varying quark masses and QCD scale will have a
substantial effect on all nuclear reaction rates as well as the binding energy
of the deuteron.

We find that the 4He abundance can be noticeably affected if the change in the
Higgs vev is a few parts per thousand or more. If ${\delta v}/{v}$ is between
$-0.03$ and $-0.01$, then the prediction fits well within the EMPRESS
calculation error bars. If one instead uses the older PDG results, then one
must (at $2\sigma$) have ${\delta v}/{v}$ less than $0.004$ and greater than
$-0.01$.

The deuterium abundance is also affected. Here, a negative value of ${\delta
v}/{v}$ makes the discrepancy with the theoretical prediction worse. A
positive value will fit within $2\sigma$ as long as it is less then $0.08$.
Together, we see that no value will be able to explain both the EMPRESS and
deuterium anomalies, but can certainly explain either one. Clearly, more
experimental results for the 4He abundance are needed as are more theoretical
studies of the deuterium abundance.

###### Acknowledgements.

We thank Albert Stadler for discussions. AKB and TMPT thank Mauro Valli for
previous collaboration on PRyMordial. The work of TMPT is supported in part by
the U.S. National Science Foundation under Grant PHY-2210283. The work of MS
is supported by the U.S. National Science Foundation under Grant PHY-2112460.
The work of AKB was supported in part by Grant No. NSF PHY-1748958 to the
Kavli Institute for Theoretical Physics (KITP), the Heising-Simons Foundation,
and the Simons Foundation (216179, LB). VK acknowledges financial support from
the Science Foundation Ireland Grant 21/PATHS/9475 (MOREHIGGS) under the SFI-
IRC Pathway Programme.

##  Appendix A The evolution of the strong coupling constant

A crude approximation to the scaling is obtained by integrating the one-loop
renormalization group equation for the strong coupling constant from a scale
well above the top quark mass down to somewhat below the charm quark mass and
then using the standard formula to deduce $\Lambda_{QCD}$. Using a mass-
independent renormalization scheme at one-loop, for $Q$ much larger than
$2m_{t}$, one has

$\frac{1}{\alpha(1\ \rm
GeV)^{2}}-\frac{1}{\alpha(Q^{2})}=\frac{1}{12\pi}(21\log 4m_{t}/Q^{2}+23\log
m_{b}^{2}/m^{2}_{t}+25\log m^{2}_{c}/m^{2}_{b}+27\log 1\ \rm
GeV^{2}/m^{2}_{c})$

Now if one multiplies all of the quark masses by $1+{\delta v}/{v}$, it is
found that

$\frac{1}{\alpha_{new}(1\ \rm GeV)^{2}}-\frac{1}{\alpha_{old}(1\ \rm
GeV)^{2}}=\frac{1}{12\pi}(21-27)*2\frac{\delta
v}{v}=-\frac{1}{\pi}\frac{\delta v}{v}$

and plugging into the standard one-loop formula for $\Lambda_{QCD}$,

$\Lambda^{2}\equiv\mu^{2}e^{-\frac{4\pi}{\beta_{0}\alpha_{s}\left(\mu^{2}\right)}},$
(2)

one finds that

$\Lambda_{new}=\mu\ \exp\frac{-6\pi}{27\alpha_{new}}=\mu\
\exp\left(\frac{-6\pi}{27\alpha_{old}}+\frac{6}{27}\frac{\delta
v}{v}\right)=\Lambda_{old}(1+\frac{2}{9}\frac{\delta v}{v})$ (3)

which is close to the results in Agrawal:1997gf ; Agrawal:1998xa . Thus, our
crude approximation yields $\Lambda_{QCD}$ scaling as $({\delta v}/{v})^{2/9}$

To be more precise, we integrate the 2 loop renormalization group equations.
The scale dependence of the strong coupling constant is controlled by the
$\beta$-function which can be expressed as a perturbative series:

$\displaystyle Q^{2}\frac{\partial}{\partial Q^{2}}\frac{\alpha_{s}}{4\pi}$
$\displaystyle=\beta\left(\alpha_{s}\right)=-\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\sum_{n=0}\left(\frac{\alpha_{s}}{4\pi}\right)^{n}\beta_{n}.$
(4)

The values of the first terms of the $\beta$-series are:

$\beta_{0}=11-\frac{2}{3}n_{f},$ (5)

at 1-loop level and

$\beta_{1}=102-\frac{38}{3}n_{f},$ (6)

at 2-loop level with $n_{f}$ the number of quark flavors active at the scale
$Q^{2}$.

The exact analytical solution to Eq. (4) is known only to $\beta_{0}$ order,
and thus we integrate numerically, defining $\Lambda_{QCD}$ to be the scale at
which $\alpha_{s}$ diverges. We find that $\Lambda_{QCD}$ scales as $({\delta
v}/{v})^{0.25}$, with the exponent varying from $0.245$ to $0.255$ over the
entire possible range of $\Lambda_{QCD}$. Thus, we use $0.25$ as the exponent.

## References

  * (1) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] [arXiv:1807.06209 [astro-ph.CO]].
  * (2) A. G. Riess, W. Yuan, L. M. Macri, D. Scolnic, D. Brout, S. Casertano, D. O. Jones, Y. Murakami, L. Breuval and T. G. Brink, et al. “A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team,” Astrophys. J. Lett. 934 (2022) no.1, L7 [arXiv:2112.04510 [astro-ph.CO]].
  * (3) M. Dainotti, B. De Simone, G. Montani, T. Schiavone and G. Lambiase, “The Hubble constant tension: current status and future perspectives through new cosmological probes,” [arXiv:2301.10572 [astro-ph.CO]].
  * (4) A. Matsumoto, M. Ouchi, K. Nakajima, M. Kawasaki, K. Murai, K. Motohara, Y. Harikane, Y. Ono, K. Kushibiki and S. Koyama, et al. “EMPRESS. VIII. A New Determination of Primordial He Abundance with Extremely Metal-poor Galaxies: A Suggestion of the Lepton Asymmetry and Implications for the Hubble Tension,” Astrophys. J. 941 (2022) no.2, 167 [arXiv:2203.09617 [astro-ph.CO]].
  * (5) O. A. Kurichin, P. A. Kislitsyn, V. V. Klimenko, S. A. Balashev and A. V. Ivanchik, “A new determination of the primordial helium abundance using the analyses of H II region spectra from SDSS,” Mon. Not. Roy. Astron. Soc. 502 (2021) no.2, 3045-3056 [arXiv:2101.09127 [astro-ph.CO]].
  * (6) T. Hsyu, R. J. Cooke, J. X. Prochaska and M. Bolte, “The PHLEK Survey: A New Determination of the Primordial Helium Abundance,” Astrophys. J. 896 (2020) no.1, 77 [arXiv:2005.12290 [astro-ph.GA]].
  * (7) M. Kawasaki and K. Murai, “Lepton asymmetric universe,” JCAP 08 (2022) no.08, 041 [arXiv:2203.09713 [hep-ph]].
  * (8) A. K. Burns, T. M. P. Tait and M. Valli, “Indications for a Nonzero Lepton Asymmetry from Extremely Metal-Poor Galaxies,” Phys. Rev. Lett. 130 (2023) no.13, 131001 [arXiv:2206.00693 [hep-ph]].
  * (9) D. Borah and A. Dasgupta, “Large neutrino asymmetry from TeV scale leptogenesis,” Phys. Rev. D 108 (2023) no.3, 035015 [arXiv:2206.14722 [hep-ph]].
  * (10) M. Escudero, A. Ibarra and V. Maura, Phys. Rev. D 107 (2023) no.3, 035024 [arXiv:2208.03201 [hep-ph]].
  * (11) P.A.M. Dirac, Nature 139 (1937) 323
  * (12) L. Bergstrom, S. Iguri and H. Rubinstein, “Constraints on the variation of the fine structure constant from big bang nucleosynthesis,” Phys. Rev. D 60 (1999), 045005 [arXiv:astro-ph/9902157 [astro-ph]].
  * (13) K. M. Nollett and R. E. Lopez, “Primordial nucleosynthesis with a varying fine structure constant: An Improved estimate,” Phys. Rev. D 66 (2002), 063507 [arXiv:astro-ph/0204325 [astro-ph]].
  * (14) M. E. Mosquera, C. G. Scoccola, S. J. Landau and H. Vucetich, “Time variation of the fine structure constant in the early universe and the Bekenstein model,” Astron. Astrophys. 478 (2008), 675 [arXiv:0707.0661 [astro-ph]].
  * (15) J. D. Barrow and J. Magueijo, “Cosmological constraints on a dynamical electron mass,” Phys. Rev. D 72 (2005), 043521 [arXiv:astro-ph/0503222 [astro-ph]].
  * (16) C. G. Scoccola, M. E. Mosquera, S. J. Landau and H. Vucetich, “Time variation of the electron mass in the early universe and the Barrow-Magueijo model,” Astrophys. J. 681 (2008), 737 [arXiv:0803.0247 [astro-ph]].
  * (17) J. Zhang and J. A. Frieman, ‘Mirror dark sector solution of the Hubble tension with time-varying fine-structure constant,” Phys. Rev. D 107 (2023) no.4, 043529 [arXiv:2211.03236 [astro-ph.CO]].
  * (18) O. Seto and Y. Toda, “Big bang nucleosynthesis constraints on varying electron mass solution to the Hubble tension,” Phys. Rev. D 107 (2023) no.8, 083512 [arXiv:2206.13209 [astro-ph.CO]].
  * (19) R. Solomon, G. Agarwal and D. Stojkovic, “Environment dependent electron mass and the Hubble constant tension,” Phys. Rev. D 105 (2022) no.10, 103536 [arXiv:2201.03127 [hep-ph]].
  * (20) L. Hart and J. Chluba, “New constraints on time-dependent variations of fundamental constants using Planck data,” Mon. Not. Roy. Astron. Soc. 474 (2018) no.2, 1850-1861 [arXiv:1705.03925 [astro-ph.CO]].
  * (21) L. Hart and J. Chluba, “Updated fundamental constant constraints from Planck 2018 data and possible relations to the Hubble tension,” Mon. Not. Roy. Astron. Soc. 493 (2020) no.3, 3255-3263 [arXiv:1912.03986 [astro-ph.CO]].
  * (22) L. Hart and J. Chluba, “Varying fundamental constants principal component analysis: additional hints about the Hubble tension,” Mon. Not. Roy. Astron. Soc. 510 (2022) no.2, 2206-2227 [arXiv:2107.12465 [astro-ph.CO]].
  * (23) O. Seto, T. Takahashi and Y. Toda, “Variation of the fine structure constant in light of recent helium abundance measurement,” Phys. Rev. D 108 (2023) no.2, 023525 [arXiv:2305.16624 [astro-ph.CO]].
  * (24) V. V. Dixit and M. Sher, “Variation of the Fermi Constant and Primordial Nucleosynthesis,” Phys. Rev. D 37 (1988), 1097
  * (25) J. Yoo and R. J. Scherrer, “Big bang nucleosynthesis and cosmic microwave background constraints on the time variation of the Higgs vacuum expectation value,” Phys. Rev. D 67 (2003), 043517 [arXiv:astro-ph/0211545 [astro-ph]].
  * (26) N. Chamoun, S. J. Landau, M. E. Mosquera and H. Vucetich, “Helium and deuterium abundances as a test for the time variation of the fine structure constant and the higgs vacuum expectation value,” J. Phys. G 34 (2007), 163-176 [arXiv:astro-ph/0508378 [astro-ph]].
  * (27) J. M. Gassner and H. Lesch, “Constraining the Spatial and Time Variation of the Higgs Vacuum Expectation Value by the Primordial He-4 Abundance,” Int. J. Theor. Phys. 47 (2008), 438-445 [arXiv:astro-ph/0606459 [astro-ph]].
  * (28) S. J. Landau, M. E. Mosquera, C. G. Scoccola and H. Vucetich, “Early Universe Constraints on Time Variation of Fundamental Constants,” Phys. Rev. D 78 (2008), 083527 [arXiv:0809.2033 [astro-ph]].
  * (29) C. M. Muller, G. Schafer and C. Wetterich, “Nucleosynthesis and the variation of fundamental couplings,” Phys. Rev. D 70, 083504 (2004) [arXiv:astro-ph/0405373 [astro-ph]].
  * (30) A. Coc, N. J. Nunes, K. A. Olive, J. P. Uzan and E. Vangioni, “Coupled Variations of Fundamental Couplings and Primordial Nucleosynthesis,” Phys. Rev. D 76, 023511 (2007) [arXiv:astro-ph/0610733 [astro-ph]].
  * (31) J. P. Kneller and G. C. McLaughlin, “BBN and Lambda(QCD),” Phys. Rev. D 68 (2003), 103508 [arXiv:nucl-th/0305017 [nucl-th]].
  * (32) V. F. Dmitriev, V. V. Flambaum and J. K. Webb, “Cosmological variation of deuteron binding energy, strong interaction and quark masses from big bang nucleosynthesis,” Phys. Rev. D 69 (2004), 063506 [arXiv:astro-ph/0310892 [astro-ph]].
  * (33) P. F. Bedaque, T. Luu and L. Platter, “Quark mass variation constraints from Big Bang nucleosynthesis,” Phys. Rev. C 83 (2011), 045803 [arXiv:1012.3840 [nucl-th]].
  * (34) G. Laverda and J. Rubio, “The Rise and Fall of the Standard-Model Higgs: Electroweak Vacuum Stability during Kination,” [arXiv:2402.06000 [hep-ph]].
  * (35) V. Agrawal, S. M. Barr, J. F. Donoghue and D. Seckel, “Viable range of the mass scale of the standard model,” Phys. Rev. D 57 (1998), 5480-5492 [arXiv:hep-ph/9707380 [hep-ph]].
  * (36) V. Agrawal, S. M. Barr, J. F. Donoghue and D. Seckel, “Anthropic considerations in multiple domain theories and the scale of electroweak symmetry breaking,” Phys. Rev. Lett. 80 (1998), 1822-1825 [arXiv:hep-ph/9801253 [hep-ph]].
  * (37) T. E. Jeltema and M. Sher, “The Triple alpha process and the anthropically allowed values of the weak scale,” Phys. Rev. D 61 (2000), 017301 [arXiv:hep-ph/9905494 [hep-ph]].
  * (38) A. K. Burns, T. M. P. Tait and M. Valli, “PRyMordial: The First Three Minutes, Within and Beyond the Standard Model,” [arXiv:2307.07061 [hep-ph]].
  * (39) S. Weinberg, “Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity,” John Wiley and Sons, 1972, ISBN 978-0-471-92567-5, 978-0-471-92567-5
  * (40) U. G. Meissner, “Low-Energy Hadron Physics from Effective Chiral Lagrangians with Vector Mesons,” Phys. Rept. 161 (1988), 213
  * (41) W. Lin and B. D. Serot, “Mid range nucleon nucleon interaction in the linear sigma model,” Nucl. Phys. A 512 (1990), 637-668
  * (42) R. L. Workman et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2022, 083C01 (2022)
  * (43) U. G. Meißner, B. C. Metsch and H. Meyer, “The electromagnetic fine-structure constant in primordial nucleosynthesis revisited,” Eur. Phys. J. A 59, no.10, 223 (2023) [arXiv:2305.15849 [hep-th]].
  * (44) N. Schöneberg, “The 2024 BBN baryon abundance update,” [arXiv:2401.15054 [astro-ph.CO]].
  * (45) A. Giri and R. J. Scherrer, “Big bang nucleosynthesis with rapidly varying G,” [arXiv:2312.06525 [astro-ph.CO]].
  * (46) H. Hong, U. Min, M. Son and T. You, “A Cosmic Window on the Dark Axion Portal,” [arXiv:2310.19544 [hep-ph]].
  * (47) F. Bianchini, G. G. di Cortona and M. Valli, “The QCD Axion: Some Like It Hot,” [arXiv:2310.08169 [hep-ph]].
  * (48) T. Chowdhury and S. Ipek, “Neutron Lifetime Anomaly and Big Bang Nucleosynthesis,” [arXiv:2210.12031 [hep-ph]].