# Optical Indistinguishability via Twinning Fields

Gerard McCaul<EMAIL_ADDRESS>Tulane University, New Orleans, LA 70118,
USA  Alexander F. King  Tulane University, New Orleans, LA 70118, USA  Denys
I. Bondar  Tulane University, New Orleans, LA 70118, USA

###### Abstract

Here we introduce the concept of the twinning field – a driving
electromagnetic pulse that induces an identical optical response from two
distinct materials. We show that for a large class of pairs of generic many-
body systems, a twinning field which renders the systems _optically
indistinguishable_ exists. The conditions under which this field exists are
derived, and this analysis is supplemented by numerical calculations of
twinning fields for both the 1D Fermi-Hubbard model, and tight-binding models
of graphene and hexagonal Boron Nitride. The existence of twinning fields may
lead to new research directions in non-linear optics, materials science, and
quantum technologies.

††preprint: APS/123-QED

#### Introduction:-

As our understanding of the physical world has progressed to mastery over it,
it has become apparent that the qualities which define a material at
equilibrium may be modified under driving. This phenomenon underpins both
quantum simulation [1] and Floquet engineering [2, 3, 4]. One of the principal
goals of quantum control theory [5] is the specification of the driving fields
necessary either to steer a system to some desired state [6, 7, 8, 9, 10, 11],
or fulfill a pre-specified condition on its expectations [12, 13].

Using tracking quantum control [14, 15, 16, 17, 18, 19, 20, 21], recent work
has demonstrated that almost arbitrary control over the optical response of a
large class of solid-state systems can be achieved [22, 2]. One consequence of
this is that two specially tailored driving fields will induce an identical
response from two distinct systems. Given the essential malleability of
quantum systems under driving, one might ask whether it is possible to fulfil
the stronger condition of obtaining identical responses using the _same_
driving field on each system. Put differently, do there exist fields for which
a pair of systems’ response are indistinguishable?

Consider two distinct systems $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$, with an
identical control field impinging on each of them (see Fig. 1). Each system
will generate an optical response $J^{(k)}(t)$, and if $J^{(1)}(t)=J^{(2)}(t)$
for all times $t$, then the driving field is what we term a ‘twinning field’,
and the systems are _optically indistinguishable_.

In the regime of linear response, this may initially appear trivial, as many
systems possess extremely similar absorption and emission spectra over a broad
range of frequencies (e.g., large organic molecules [24]). Indeed, such is the
closeness of these systems’ response that quantum control [25, 26, 27, 28, 29,
30, 31] (including tracking control [32]) must be exploited to accurately
detect these systems. Of course, similar is not _identical_ , and purely
linear response would require identical susceptibilities for identical
responses. In general however, materials also have a non-linear component to
differentiate them (see, e.g., [33]). In fact many important phenomena – e.g.
high harmonic generation [34, 35, 36], the workhorse of attosecond physics
[37, 38, 39] – explicitly rely on optical non-linearities [40].

Figure 1: Ordinarily, distinct systems will have different responses to the
same driving field. A twinning field relates a pair of systems as the field
under which the optical response $J^{(k)}(t)$ of each system is identical.

For this reason, true optical indistinguishability must be considered in the
context of the non-linear response that arises in a fully quantum treatment of
materials. Some preliminary hints of this indistinguishability have been
observed experimentally, specifically in the non-uniqueness of the
parametrisation of second-order nonlinear spectra [41]. To date however, a
theoretical justification for such results has been lacking. In this
manuscript we address this issue and present a framework for achieving driven
indistinguishability. The main result is a demonstration that for any pair of
generic many-electron systems on a lattice, there always exists a twinning
field which will elicit an identical response from each system. Furthermore,
the conditions under which this field is unique are established. This
framework is then extended by deriving a general twinning field which renders
an arbitrary expectation identical between systems. Finally we discuss the
physical implications of twinning fields, and their potential utility.

#### Results:-

Here we outline the derivation of the twinning field for 1D systems, while the
more general case is discussed in the supplementary material [42] which
includes Refs. [4, 5, 3]. Here we consider two many-electron systems on a
lattice, labeled by $k=1,2$. Each has a potential $\hat{U}^{(k)}$ due to
electron-electron interactions. Both systems are excited by an identical laser
pulse, described under the dipole approximation by the Peierls phase $\Phi(t)$
[46, 47]. Such systems’ evolution will then be determined by the Hamiltonian
(in atomic units) [2, 48]:

$\displaystyle\hat{H}^{(k)}(t)=-t^{(k)}_{0}\sum_{j,\sigma}\text{$\left({\rm
e}^{-i\Phi\left(t\right)}\hat{c}_{j\sigma}^{\dagger}\hat{c}_{j+1\sigma}+\rm{h.c.}\right)$}+\hat{U}^{(k)},$
(1)

where $\hat{c}_{j\sigma}$ is the fermionic annhilation operator (acting on the
appropriate system) for site $j$ and spin $\sigma$, satisfying the
anticommutation relation
$\\{\hat{c}^{\dagger}_{j\sigma},\hat{c}_{j^{\prime}\sigma^{\prime}}\\}=\delta_{\sigma\sigma^{\prime}}\delta_{jj^{\prime}}$,
while $t^{(k)}_{0}$ is the hopping parameter describing the kinetic energy of
the electrons.

Consider a typical example of an optically driven current. The current
operator $\hat{J}^{(k)}$ is defined from a continuity equation for the
electron density [22, 2]. Provided all number operators
$\hat{n}_{j\sigma}=\hat{c}^{\dagger}_{j\sigma}\hat{c}_{j\sigma}$ commute with
$\hat{U}^{(k)}$, each system’s current operator has the form [49]:

$\hat{J}^{(k)}(t)=-iat^{(k)}_{0}\sum_{j,\sigma}\left({\rm
e}^{-i\Phi\left(t\right)}\hat{c}_{j\sigma}^{\dagger}\hat{c}_{j+1\sigma}-{\rm
h.c.}\right),$ (2)

where $a^{(k)}$ is the lattice constant. It is important to note that the
current expectation
$J^{(k)}(t)=\braket{\psi_{k}(t)}{\hat{J}^{(k)}(t)}{\psi_{k}(t)}$ depends only
implicitly on $\hat{U}^{(k)}$ through the evolution of $\ket{\psi_{k}(t)}$,
significantly simplifying expressions.

Having dispensed with this preamble, we come to our main topic of
investigation. For two systems $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$ with
potentials $\hat{U}^{(1)}$ and $\hat{U}^{(2)}$, does there exist a twinning
field $\Phi_{\tau}(t)$ such that $J^{(1)}(t)=J^{(2)}(t)$, making the response
of one system indistinguishable from the other?

To establish the existence of this field, we first express the nearest-
neighbour expectation of each system in a polar form

$\displaystyle\hat{K}$
$\displaystyle=\sum_{j,\sigma}\hat{c}_{j\sigma}^{\dagger}\hat{c}_{j+1\sigma},$
(3) $\displaystyle K(\psi_{k})$
$\displaystyle=\left\langle\psi_{k}(t)\left|\hat{K}\right|\psi_{k}(t)\right\rangle=R\left(\psi_{k}\right){\rm
e}^{i\theta\left(\psi_{k}\right)}.$ (4)

Note that in both this and later expressions, the argument $\psi_{k}$
indicates that the expression is a functional of
$\ket{\psi_{k}}\equiv\ket{\psi_{k}(t)}$. We emphasise that this functional
will have a well defined value for any state that has been obtained through
evolution under the Hamiltonian given in Eq. (1). Using this, we may express
the response expectation $J^{(k)}$ directly as

$\displaystyle J^{(k)}\left(t\right)=$ $\displaystyle-
iat_{0}R\left(\psi_{k}\right)\left({\rm
e}^{-i\left[\Phi\left(t\right)-\theta\left(\psi_{k}\right)\right]}-{\rm
e}^{i\left[\Phi\left(t\right)-\theta\left(\psi_{k}\right)\right]}\right)$
$\displaystyle=$
$\displaystyle-2at_{0}R\left(\psi_{k}\right)\sin(\Phi\left(t\right)-\theta\left(\psi_{k}\right)).$
(5)

It is straightforward to equate the currents $J^{(1)}(t)$ and $J^{(2)}(t)$ and
obtain an expression for the twinning field $\Phi_{\tau}$ in terms of the
expectations of the two systems:

$\displaystyle\Phi_{\tau}(t)$
$\displaystyle=\arctan\left(\xi(\psi_{1},\psi_{2})\right)$ (6)
$\displaystyle\xi(\psi_{1},\psi_{2})$ $\displaystyle=\frac{\lambda
R(\psi_{1})\sin(\theta(\psi_{1}))-R(\psi_{2})\sin(\theta(\psi_{2}))}{\lambda
R(\psi_{1})\cos(\theta(\psi_{1}))-R(\psi_{2})\cos(\theta(\psi_{2}))},$ (7)

where $\lambda=\frac{a^{(1)}t^{(1)}_{0}}{a^{(2)}t^{(2)}_{0}}$. Critically, in
this 1D case one is able to obtain a closed form for $\Phi_{\tau}(t)$, such
that the existence of this field can be assessed purely by considering its
right hand side. Given both the range and domain of $\arctan$ extends over the
reals and $\xi(\psi_{1},\psi_{2})$ is real by definition (and has a definite
value for any pair of states), we can immediately conclude that _a twinning
field between any two systems described by Eq. ( 1) always exists._

An important caveat to this statement is that the predicted twinning field may
be identically zero depending on the initial states of the twinned systems.
For example, if we attempt to twin two systems of non-interacting electrons
($\hat{U}^{(k)}=0$), then $\hat{K}$ commutes with the Hamiltonian and
$K(\psi_{k})$ is constant. If the systems are evolved from their ground state,
$\theta(\psi_{k})=0$, and by Eq. (6), $\Phi_{\tau}(t)=0$. This scenario is
consistent with the impossibility of twinning fields in linear optics, and can
be avoided by having at least one of the system pair have a non-zero potential
and hence a nonlinear response. Furthermore, while an equation for
$\Phi_{\tau}$ can still be obtained in higher dimensions, in general it will
not be of a closed form, and therefore a twinning field is not guaranteed to
exist. The additional requirements for a twinning field to exist in this
scenario are detailed in the supplementary material [42].

#### Illustrations:-

Here we provide examples of twinning fields for systems described by the
Fermi-Hubbard model of strongly interacting electrons. In this case, each
system has an onsite potential described by [48, 49]:

$\hat{U}^{(k)}=U^{(k)}\sum_{j}\hat{c}_{j\uparrow}^{\dagger}\hat{c}_{j\uparrow}\hat{c}_{j\downarrow}^{\dagger}\hat{c}_{j\downarrow}$
(8)

where $U^{(k)}$ parametrises the energy of the electron-electron repulsion.
Systems’ equilibrium properties are determined by the ratio
$U^{(k)}/t^{(k)}_{0}$. Despite the simplicity of the potential, this model is
rich in nontrivial behaviour, including topological [50, 51] and
superconducting phases in 2D [52, 53]. The Fermi-Hubbard model is
computationally challenging and a complete understanding of its dynamics is
believed to require a quantum computer [54]. It also exhibits a highly non-
linear optical response [55, 49, 56], and therefore provides a suitable
platform for numerical calculations of twinning fields.

Here we consider an $L=10$ site chain with periodic boundary conditions and an
average of one electron per site. For the sake of simplicity, in both systems
we use the lattice constants $a^{(1)}=a^{(2)}=4$ Å, with a hopping parameter
of $t_{0}\equiv t^{(1)}_{0}=t_{0}^{(2)}=0.52$ eV. To avoid the trivial
solution of $\Phi_{\tau}(t)=0$, each system (initially in the ground state) is
first pumped by a single cycle of a transform-limited field. Specifically, an
enveloped sine-wave is used with an amplitude of $E_{0}=10$ MV/cm and
frequency $\omega_{0}=32.9$ THz. All calculations were performed using exact
diagonalization via the QuSpin Python package [57, 58].

Figure 2: Top panel: twinning fields and the accompanying optical response for
three pairs of systems, following an initial pump pulse. Bottom panel: The
current resulting from the application of the twinning field, which is
identical in each pair of systems. In all cases $U^{(1)}=0.5t_{0}$, while
$U^{(2)}=U^{(1)}+\Delta$, where $\Delta$ is varied for each of the three
pairs. ${\rm a^{\prime}.u.}$ are atomic units with energy normalised to
$t_{0}$.

Figure 2 shows examples of calculated twinning fields and the accompanying
responses they generate for several pairs of systems. These pairs are
parametrised by $\Delta$, using $U^{(1)}=0.5t_{0}$, while
$U^{(2)}=U^{(1)}+\Delta$. Applying the twinning field calculated at each time,
we find the current in each pair of systems is identical, as expected.

Figure 3: Both graphene and hBN can be modelled with the tight-binding
approximation set on a hexagonal bipartite lattice, defined by the positions
of the $B$ sublattice atoms relative to those on the $A$ sublattice. These are
characterised by the three nearest neighbour vectors $\mathbf{a}_{j}$, with
each being having an angular separation of $120^{\circ}$, and length 2.5$\AA$

A more concretely physical example is to twin two commonly studied materials -
graphene and hexagonal Boron Nitride (hBN). Both of these structures have a
bipartite lattice structure, as shown in Fig.3. Critically, that lattice
constant for both systems is almost identical, and can hence be modelled with
the identical value
$\left|\mathbf{a}^{(gr)}_{j}\right|=\left|\mathbf{a}^{(hBN)}_{j}\right|=2.5\AA$
[59, 60]. This is critical for the existence of a twinning field, as in higher
dimensions the ability to twin a pair of systems can only be guaranteed when
they share the same lattice structure.

Both materials are well modelled by the tight-binding approximation [61], with
an onsite potential of the same form as Eq. (8).
$t^{(gr)}_{0}=t^{(hBN)}_{0}=2.7$eV, and for the graphene carbon atoms
$U_{c}=0$, while for hBN $U_{B}=3.3$eV and $U_{N}=-1.4$eV for Boron and
Nitrogen atoms respectively [62, 63]. While it is possible to also include
next-to-nearest hopping, the relative strength of this compared to nearest-
neighbour hopping is only $\sim 5$% [64, 65], and we therefore neglect it for
calculational simplicity. Further information on both the precise Hamiltonian
describing these systems, and the derived twinning field equations may be
found in the supplementary information [42].

Simulations are performed using $L=12$ sites with periodic boundaries, and the
systems are again prepared via the application of a pumping field. In this
case, the polarisation of this initial pump is of great consequence, and in
order to generate a physically realisable twinning field, it must be aligned
with one of the nearest neighbour vectors. Fig. 4 shows an example of this,
with the initial pump pulse and subsequent twinning field aligned along the
$\mathbf{a}_{1}$ direction.

Figure 4: Top panel: Twinning field generated after a pump pulse in the
$\mathbf{a}_{1}$/$x$ direction. The symmetry of the lattice guarantees that
the resultant twinning field is _also_ purely in the $x$ direction. Bottom
panel: the overall current $J(t)$ in the $x$ direction due to application of
the twinning field. As expected, this is identical in both systems.

#### Discussion:-

We have introduced a non-linear optical phenomenon where some pairs of quantum
systems have a twinning field which generates an identical response in each
individual system. In 1D the necessary conditions for a twinning field to
exist are rather general, but in higher dimensions the two systems must
possess a high degree of similarity in their lattice structures for the
existence of a twinning field to be guaranteed. Conditions for the uniqueness
of this field were derived, and numerical calculations provided examples of
twinning fields in a Fermi-Hubbard system.

It is instructive to compare optical indistinguishability with the anti-
haecceitism [66, 67] of quantum particles. The latter is responsible for both
the Fermi-Dirac and Bose-Einstein distributions [68, 69] (as well as the
resolution of Gibbs’ paradox [70, 71]), and is an intrinsic and immutable
property of said particles. For this reason it has been commonly assumed that
systems governed by distinct Hamiltonians will be distinguishable from each
other. Indeed, the effectiveness of spectroscopy is predicated on the notion
that a material can be uniquely identified from its spectral response [72].
The existence of twinning fields gives the lie to this assumption however,
demonstrating that indistinguishability can arise as an emergent property
under driving.

It may be tempting to think of twinning fields purely as an act of deception,
where an unscrupulous salesman could use the technique to pass off a cheap and
nasty material as something more costly. In fact, the analysis presented here
demonstrates that it is a trick requiring highly specific conditions to be
repeated, and such a fraud can be defeated by an arbitrary modification to the
example driving field. Indeed, even if the twinning field is non-unique, the
field up to the point that Lipschitz continuity (see [42]) is violated _will_
be unique. Any field that is distinct from this initial trajectory will
therefore be _guaranteed to produce a response distinguishing the two
systems_. This has the important consequence of ensuring that techniques
designed to discriminate between similar systems are well founded [27, 28,
32].

Of course, the existence of twinning fields forces one to consider both their
feasibility and wider utility. While the twinning fields calculated here
appear to be rather broadband pulses, the rapid improvement in both intensity
and bandwidth of laboratory laser sources [73, 74] – combined with the fact
that similarly tailored tracking control fields can be well approximated by a
few distinct frequencies [2] – suggests twinning fields may be experimentally
realisable with current technology. In fact, deep learning networks have
recently been employed to experimentally determine the driving field required
to generate a desired response in a material, in a manner that is robust to
noise [75]. It is likely that such techniques could be similarly applied for
the practical calculation of twinning fields.

One potential application of these fields is to characterise the effect of
interactions between systems. Applying the twinning field calculated for a
non-interacting pair, it would be possible to identify the additional current
generated by each system due to its interaction with the other. Twinning
fields may also provide a method for creating alternative realisations of
metamaterials [76] when a specific response is required. Given a field and a
metamaterial’s nearest-neighbour expectation, Eq. (6) can be used to calculate
the properties of a different material producing the same effect. Consequently
the search for cheaper alternative components in quantum technologies may be
aided by twinning fields.

Naturally, there remain a number of unanswered questions. For instance, a
given twinning field relates a pair of systems, but is this pair unique? Put
differently, are there triplets, or $n$-tuplets of systems which exhibit
optical indistinguishability? In the case of a non-unique twinning field, what
are the physical consequences of choosing one solution over an another? Such
questions may merit further investigation, as understanding these secondary
properties provides both challenges and opportunities to illuminate the
principles upon which driven systems operate.

_Acknowledgments:_ The authors are grateful to the anonymous referee for
comments which have strengthened the paper. The authors would also like to
thank Alicia Magann their helpful comments on the manuscript. This work has
been generously supported by Army Research Office (ARO) (grant
W911NF-19-1-0377; program manager Dr. James Joseph). The views and conclusions
contained in this document are those of the authors and should not be
interpreted as representing the official policies, either expressed or
implied, of ARO or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein. This work has also been partially supported by by
the W. M. Keck Foundation.

## References

  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * Oka and Kitamura [2019] T. Oka and S. Kitamura, Floquet engineering of quantum materials, Annual Review of Condensed Matter Physics 10, 387 (2019), https://doi.org/10.1146/annurev-conmatphys-031218-013423 .
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics 64, 139 (2015), https://doi.org/10.1080/00018732.2015.1055918 .
  * Claeys  _et al._ [2019] P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov, Floquet-engineering counterdiabatic protocols in quantum many-body systems, Phys. Rev. Lett. 123, 090602 (2019).
  * Werschnik and Gross [2007] J. Werschnik and E. K. Gross, Quantum optimal control theory, J. Phys. B 40 (2007).
  * Serban  _et al._ [2005] I. Serban, J. Werschnik, and E. K. Gross, Optimal control of time-dependent targets, Phys. Rev. A 71, 053810 (2005).
  * Doria  _et al._ [2011] P. Doria, T. Calarco, and S. Montangero, Optimal control technique for many-body quantum dynamics, Phys. Rev. Lett. 106, 190501 (2011).
  * Peirce  _et al._ [1988] A. P. Peirce, M. A. Dahleh, and H. Rabitz, Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications, Phys. Rev. A 37, 4950 (1988).
  * Winterfeldt  _et al._ [2008] C. Winterfeldt, C. Spielmann, and G. Gerber, Colloquium: Optimal control of high-harmonic generation, Rev. Mod. Phys. 80, 117 (2008).
  * Glaser  _et al._ [2015] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Training schrödinger’s cat: quantum optimal control, Eur. Phys. J. D 69, 279 (2015).
  * Magann  _et al._ [2021] A. B. Magann, C. Arenz, M. D. Grace, T.-S. Ho, R. L. Kosut, J. R. McClean, H. A. Rabitz, and M. Sarovar, From pulses to circuits and back again: A quantum optimal control perspective on variational quantum algorithms, PRX Quantum 2, 010101 (2021).
  * Kosloff  _et al._ [1992] R. Kosloff, A. D. Hammerich, and D. Tannor, Excitation without demolition: Radiative excitation of ground-surface vibration by impulsive stimulated Raman scattering with damage control, Phys. Rev. Lett. 69, 2172 (1992).
  * Bartana  _et al._ [1993] A. Bartana, R. Kosloff, and D. J. Tannor, Laser cooling of molecular internal degrees of freedom by a series of shaped pulses, J. Chem. Phys. 99, 196 (1993).
  * Ong  _et al._ [1984] C. K. Ong, G. M. Huang, T. J. Tarn, and J. W. Clark, Invertibility of quantum-mechanical control systems, Mathematical systems theory 17, 335 (1984).
  * Clark  _et al._ [1985] J. W. Clark, C. K. Ong, T. J. Tarn, and G. M. Huang, Quantum nondemolition filters, Mathematical systems theory 18, 33 (1985).
  * Rothman  _et al._ [2005] A. Rothman, T.-S. Ho, and H. Rabitz, Observable-preserving control of quantum dynamics over a family of related systems, Phys. Rev. A 72, 023416 (2005).
  * Magann  _et al._ [2018] A. Magann, T.-S. Ho, and H. Rabitz, Singularity-free quantum tracking control of molecular rotor orientation, Phys. Rev. A 98, 043429 (2018).
  * Caneva  _et al._ [2011] T. Caneva, T. Calarco, and S. Montangero, Chopped random-basis quantum optimization, Phys. Rev. A 84, 022326 (2011).
  * Campos  _et al._ [2017] A. G. Campos, D. I. Bondar, R. Cabrera, and H. A. Rabitz, How to Make Distinct Dynamical Systems Appear Spectrally Identical, Phys. Rev. Lett. 118, 083201 (2017).
  * Zhu and Rabitz [2003] W. Zhu and H. Rabitz, Quantum control design via adaptive tracking, J. Chem. Phys. 119, 3619 (2003).
  * Zhu  _et al._ [1999] W. Zhu, M. Smit, and H. Rabitz, Managing singular behavior in the tracking control of quantum dynamical observables, J. Chem. Phys. 110, 1905 (1999).
  * McCaul  _et al._ [2020a] G. McCaul, C. Orthodoxou, K. Jacobs, G. H. Booth, and D. I. Bondar, Driven imposters: Controlling expectations in many-body systems, Phys. Rev. Lett. 124, 183201 (2020a).
  * McCaul  _et al._ [2020b] G. McCaul, C. Orthodoxou, K. Jacobs, G. H. Booth, and D. I. Bondar, Controlling arbitrary observables in correlated many-body systems, Phys. Rev. A 101, 053408 (2020b).
  * Lichtman and Conchello [2005] J. W. Lichtman and J.-A. Conchello, Fluorescence microscopy, Nat. Methods 2, 910 (2005).
  * Brixner  _et al._ [2001] T. Brixner, N. Damrauer, P. Niklaus, and G. Gerber, Photoselective adaptive femtosecond quantum control in the liquid phase, Nature 414, 57 (2001).
  * Goun  _et al._ [2016] A. Goun, D. I. Bondar, A. O. Er, Z. Quine, and H. A. Rabitz, Photonic reagents for concentration measurement of fluorescent proteins with overlapping spectra, Sci. Rep. 6, 25827 (2016).
  * Li  _et al._ [2002] B. Li, G. Turinici, V. Ramakrishna, and H. Rabitz, Optimal dynamic discrimination of similar molecules through quantum learning control, J. Phys. Chem. B 106, 8125 (2002).
  * Li  _et al._ [2005] B. Li, H. Rabitz, and J.-P. Wolf, Optimal dynamic discrimination of similar quantum systems with time series data, J. Chem. Phys 122, 154103 (2005).
  * Rondi  _et al._ [2012] A. Rondi, L. Bonacina, A. Trisorio, C. Hauri, and J.-P. Wolf, Coherent manipulation of free amino acids fluorescence, Phys. Chem. Chem. Phys. 14, 9317 (2012).
  * Roslund  _et al._ [2011] J. Roslund, M. Roth, L. Guyon, V. Boutou, F. Courvoisier, J.-P. Wolf, and H. Rabitz, Resolution of strongly competitive product channels with optimal dynamic discrimination: Application to flavins, J. Chem. Phys. 134, 034511 (2011).
  * Roth  _et al._ [2009] M. Roth, L. Guyon, J. Roslund, V. Boutou, F. Courvoisier, J.-P. Wolf, and H. Rabitz, Quantum control of tightly competitive product channels, Phys. Rev. Lett. 102, 253001 (2009).
  * Magann  _et al._ [2020] A. B. Magann, G. McCaul, H. A. Rabitz, and D. I. Bondar, Scalable single-shot mixture characterization (2020), arXiv:2010.13859 .
  * Drobizhev  _et al._ [2011] M. Drobizhev, N. S. Makarov, S. E. Tillo, T. E. Hughes, and A. Rebane, Two-photon absorption properties of fluorescent proteins, Nature methods 8, 393 (2011).
  * Ghimire  _et al._ [2012] S. Ghimire, A. D. Dichiara, E. Sistrunk, G. Ndabashimiye, U. B. Szafruga, A. Mohammad, P. Agostini, L. F. Dimauro, and D. A. Reis, Generation and propagation of high-order harmonics in crystals, Phys. Rev. A 85, 043836 (2012).
  * Ghimire  _et al._ [2011] S. Ghimire, A. D. Dichiara, E. Sistrunk, P. Agostini, L. F. Dimauro, and D. A. Reis, Observation of high-order harmonic generation in a bulk crystal, Nature Physics 7, 138 (2011).
  * Murakami  _et al._ [2018] Y. Murakami, M. Eckstein, and P. Werner, High-Harmonic Generation in Mott Insulators, Phys. Rev. Lett. 121, 057405 (2018).
  * Corkum and Krausz [2007] P. B. Corkum and F. Krausz, Attosecond science, Nature Physics 3, 381 (2007).
  * Krausz and Ivanov [2009] F. Krausz and M. Ivanov, Attosecond physics, Rev. Mod. Phys. 81, 163 (2009).
  * Li  _et al._ [2020] J. Li, J. Lu, A. Chew, S. Han, J. Li, Y. Wu, H. Wang, S. Ghimire, and Z. Chang, Attosecond science based on high harmonic generation from gases and solids, Nature Communications 11, 10.1038/s41467-020-16480-6 (2020).
  * Boyd [2008] R. W. Boyd,  _Nonlinear optics_ , 3rd ed. (Academic Press, Amsterdam ; Boston, 2008).
  * Busson and Tadjeddine [2009] B. Busson and A. Tadjeddine, Non-uniqueness of parameters extracted from resonant second-order nonlinear optical spectroscopies, Journal of Physical Chemistry C 113, 21895 (2009).
  * [42] See Supplemental Material at [URL will be inserted by publisher] for a full derivation of the twinning field.
  * Nagle [2011] R. K. Nagle,  _Fundamentals of Differential Equations and Boundary Value Problems (6th Edition) (Featured Titles for Differential Equations)_ (Pearson, 2011).
  * Folland [2007] G. B. Folland,  _Real Analysis: Modern Techniques and Their Applications_ (Wiley, 2007).
  * Jha  _et al._ [2009] A. Jha, V. Beltrani, C. Rosenthal, and H. Rabitz, Multiple solutions in the tracking control of quantum systems, J. Phys. Chem. A 113, 7667 (2009).
  * Peierls [1933] R. Peierls, Zur theorie des diamagnetismus von leitungselektronen, Zeitschrift für Physik 80, 763 (1933).
  * Nocera  _et al._ [2017] A. Nocera, A. Polkovnikov, and A. E. Feiguin, Unconventional fermionic pairing states in a monochromatically tilted optical lattice, Phys. Rev. A 95, 023601 (2017).
  * Essler [2005] F. H. L. Essler,  _The One-Dimensional Hubbard Model_ (Cambridge University Press, 2005).
  * Silva  _et al._ [2018] R. E. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M. Ivanov, High-harmonic spectroscopy of ultrafast many-body dynamics in strongly correlated systems, Nature Photonics 12, 266 (2018).
  * Le  _et al._ [2020] N. H. Le, A. J. Fisher, N. J. Curson, and E. Ginossar, Topological phases of a dimerized fermi–hubbard model for semiconductor nano-lattices, npj Quantum Information 6, 10.1038/s41534-020-0253-9 (2020).
  * Jia  _et al._ [2020] L. Jia, Z. Zhang, D. Z. Yang, M. S. Si, and G. P. Zhang, Probing magnetic configuration-mediated topological phases via high harmonic generation in ${\mathrm{mnbi}}_{2}{\mathrm{te}}_{4}$, Phys. Rev. B 102, 174314 (2020).
  * Maier  _et al._ [2005] T. A. Maier, M. Jarrell, T. C. Schulthess, P. R. C. Kent, and J. B. White, Systematic study of $d$-wave superconductivity in the 2d repulsive hubbard model, Phys. Rev. Lett. 95, 237001 (2005).
  * Gull  _et al._ [2013] E. Gull, O. Parcollet, and A. J. Millis, Superconductivity and the pseudogap in the two-dimensional hubbard model, Phys. Rev. Lett. 110, 216405 (2013).
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Rev. Mod. Phys. 92, 015003 (2020).
  * Keimer and Moore [2017] B. Keimer and J. E. Moore, The physics of quantum materials, Nature Physics 13, 1045 (2017).
  * Tokura and Nagaosa [2018] Y. Tokura and N. Nagaosa, Nonreciprocal responses from non-centrosymmetric quantum materials, Nature Communications 9, 10.1038/s41467-018-05759-4 (2018).
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems. Part II: bosons, fermions and higher spins, SciPost Phys. 7, 20 (2019).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems part I: spin chains, SciPost Phys. 2, 003 (2017).
  * Giovannetti  _et al._ [2007] G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations, Phys. Rev. B 76, 073103 (2007).
  * Lin and Ni [2020] X. Lin and J. Ni, Symmetry breaking in the double moiré superlattices of relaxed twisted bilayer graphene on hexagonal boron nitride, Phys. Rev. B 102, 035441 (2020).
  * Reich  _et al._ [2002] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, Tight-binding description of graphene, Phys. Rev. B 66, 035412 (2002).
  * Sławińska  _et al._ [2010] J. Sławińska, I. Zasada, and Z. Klusek, Energy gap tuning in graphene on hexagonal boron nitride bilayer system, Phys. Rev. B 81, 155433 (2010).
  * Moon and Koshino [2014] P. Moon and M. Koshino, Electronic properties of graphene/hexagonal-boron-nitride moiré superlattice, Phys. Rev. B 90, 155406 (2014).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Deacon  _et al._ [2007] R. S. Deacon, K.-C. Chuang, R. J. Nicholas, K. S. Novoselov, and A. K. Geim, Cyclotron resonance study of the electron and hole velocity in graphene monolayers, Phys. Rev. B 76, 081406 (2007).
  * Saunders [2013] S. Saunders,  _Indistinguishability_ (Oxford University Press, 2013).
  * Cowling [2016] S. Cowling, Haecceitism, in  _The Stanford Encyclopedia of Philosophy_ , edited by E. N. Zalta (Metaphysics Research Lab, Stanford University, 2016) fall 2016 ed.
  * Dirac [1981] P. A. M. Dirac,  _The principles of quantum mechanics_ (Clarendon Press, Oxford England, 1981).
  * Ziff  _et al._ [1977] R. M. Ziff, G. E. Uhlenbeck, and M. Kac, The ideal bose-einstein gas, revisited, Physics Reports 32, 169 (1977).
  * Jaynes [1965] E. T. Jaynes, Gibbs vs Boltzmann Entropies, Am. J. Phys. 33, 391 (1965).
  * Jaynes [1992] E. T. Jaynes, The Gibbs paradox, Maximum entropy and bayesian methods 17, 1 (1992), arXiv:1003.0179 .
  * Workman [1998] J. Workman,  _Applied spectroscopy : a compact reference for practitioners_ (Academic Press, San Diego, 1998).
  * Elu  _et al._ [2020] U. Elu, L. Maidment, L. Vamos, F. Tani, D. Novoa, M. H. Frosz, V. Badikov, D. Badikov, V. Petrov, P. S. J. Russell, _et al._ , Seven-octave high-brightness and carrier-envelope-phase-stable light source, Nat. Photonics , 1 (2020).
  * Lesko  _et al._ [2020] D. Lesko, H. Timmers, S. Xing, A. Kowligy, A. J. Lind, and S. A. Diddams, A six-octave optical frequency comb from a scalable few-cycle erbium fiber laser, arXiv preprint arXiv:2005.13673 (2020).
  * Lohani  _et al._ [2019] S. Lohani, E. M. Knutson, W. Zhang, and R. T. Glasser, Dispersion characterization and pulse prediction with machine learning, OSA Continuum 2, 3438 (2019).
  * Lapine  _et al._ [2014] M. Lapine, I. V. Shadrivov, and Y. S. Kivshar, Colloquium: Nonlinear metamaterials, Rev. Mod. Phys. 86, 1093 (2014).

Supplemental Material:Optical Indistinguishability via Twinning Fields

##  I Twinning Field Derivation

###  I.1 General Hamiltonian

Here we provide a full derivation for the twinning field, beginning with the
continuum, non-relativistic Hamiltonian for a driven electronic system:

$\displaystyle\hat{\mathcal{H}}=\sum_{\sigma}\int\frac{{\rm
d^{3}}x}{2}\hat{\psi}_{\sigma}^{\dagger}(\mathbf{r})\left[-\nabla^{2}+V(\mathbf{r})+\mathbf{r}\cdot\mathbf{E}(\mathbf{r},t)\right]\hat{\psi}_{\sigma}(\mathbf{r})+\sum_{\sigma\sigma^{\prime}}\int\frac{{\rm
d}^{3}x{\rm
d}^{3}x^{\prime}}{2}\hat{\psi}_{\sigma^{\prime}}^{\dagger}(\mathbf{r}^{\prime})\hat{\psi}_{\sigma}^{\dagger}(\mathbf{r})U_{\sigma\sigma^{\prime}}(\mathbf{r},\mathbf{r}^{\prime})\hat{\psi}_{\sigma}(\mathbf{r})\hat{\psi}_{\sigma^{\prime}}(\mathbf{r}^{\prime})$
(S1)

Here the $\hat{\psi}_{\sigma}(\mathbf{r})$ are standard fermionic field
operators satisfying
$\left\\{\hat{\psi}_{\sigma^{\prime}}^{\dagger}(\mathbf{r}^{\prime}),\hat{\psi}_{\sigma}(\mathbf{r})\right\\}=\delta_{\sigma\sigma^{\prime}}\delta(\mathbf{r}-\mathbf{r}^{\prime})$.
$V(\mathbf{r})$ is an external potential, while
$U_{\sigma\sigma^{\prime}}(\mathbf{r},\mathbf{r}^{\prime})$ is the two-body
interaction potential, and the system is coupled to the driving electric field
via the dipole approximation $\mathbf{r}\cdot\mathbf{E}(\mathbf{r},t)$. Note
that in principle it is possible to include higher order interactions without
materially affecting the derivation, but in the interests of clarity and
concision, we do not include these terms. For similar reasons, we specify that
the driving field be homogeneous, i.e.
$\mathbf{E}(\mathbf{r},t)\equiv\mathbf{E}(t)$.

In order to establish a direct relationship between a given expectation and
the driving field, the system must be discretised in time or space. Performing
a spatial discretisation, we obtain the general form:

$\hat{H}=\sum_{\mathbf{r}\mathbf{r}^{\prime},\sigma}t(\mathbf{r},\mathbf{r}^{\prime})\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{\mathbf{r}^{\prime}\sigma}+\sum_{\mathbf{r}\mathbf{r}^{\prime},\sigma\sigma^{\prime}}U_{\sigma\sigma^{\prime}}(\mathbf{r},\mathbf{r}^{\prime})\hat{n}_{\mathbf{r}\sigma}\hat{n}_{\mathbf{r}^{\prime}\sigma^{\prime}}+\sum_{\mathbf{r},\sigma}\hat{n}_{\mathbf{r}\sigma}\mathbf{r}\cdot\mathbf{E}(t).$
(S2)

In this case one body terms have been collected into (the symmetric function)
$t(\mathbf{r},\mathbf{r}^{\prime})$, and the field operators have been
replaced by the lattice annihilation operators
$\\{\hat{c}^{\dagger}_{\mathbf{r}\sigma},\hat{c}_{\mathbf{r}^{\prime}\sigma^{\prime}}\\}=\delta_{\sigma\sigma^{\prime}}\delta_{\mathbf{r}\mathbf{r}^{\prime}}$.
Finally we have defined the number operator
$\hat{n}_{\mathbf{r}\sigma}=\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{\mathbf{r}\sigma}$.
To proceed further we perform the following time dependent unitary
transformation of the Hamiltonian:

$\displaystyle\hat{H}$
$\displaystyle\to\hat{R}^{\dagger}(t)\hat{H}\hat{R}(t)-i\hat{R}^{\dagger}(t)\frac{\partial\hat{R}(t)}{\partial
t}$ (S3) $\displaystyle\hat{R}(t)$ $\displaystyle={\rm
e}^{i\mathbf{r}\cdot\mathbf{A}(t)},\ \ \mathbf{A}(t)=\int^{t}_{-\infty}{\rm
d}t^{\prime}\ \mathbf{E}(t),$ (S4)

from which we obtain:

$\hat{H}=\sum_{\mathbf{r}\mathbf{r}^{\prime},\sigma}\left({\rm
e}^{-i(\mathbf{r}-\mathbf{r}^{\prime})\cdot\mathbf{A}(t)}t(\mathbf{r},\mathbf{r}^{\prime})\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{\mathbf{r}^{\prime}\sigma}+{\rm
h.c.}\right)+\sum_{\mathbf{r}\mathbf{r}^{\prime},\sigma\sigma^{\prime}}U_{\sigma\sigma^{\prime}}(\mathbf{r},\mathbf{r}^{\prime})\hat{n}_{\mathbf{r}\sigma}\hat{n}_{\mathbf{r}^{\prime}\sigma^{\prime}}.$
(S5)

At this point, it is helpful to specify the form of the one-body term
$t(\mathbf{r},\mathbf{r}^{\prime})$. For the sake of simplicity, we take this
to be a nearest neighbour interaction, but note it is possible to include
additional next-to-nearest neighbour (and even longer range) hopping terms
without affecting our ability to derive a closed form for the twinning field.
In the systems we consider here, such terms will be small relative to the
nearest neighbour hopping, and their presence would further clutter equations
already suffering from index overburdening.

If each site has $M$ neighbours, described by the relative vectors
$\mathbf{a}_{j}$, then we express the hopping as:

$t(\mathbf{r},\mathbf{r}^{\prime})=-t_{0}\sum_{j=1}^{j=M}\delta_{(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{a}_{j}},$
(S6)

which yields the Hamiltonian

$\hat{H}=-t_{0}\sum_{\mathbf{r},\sigma}\sum_{j}\left({\rm
e}^{-i\Phi_{j}(t)}\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{(\mathbf{r}+\mathbf{a}_{j})\sigma}+{\rm
h.c.}\right)+\hat{U},$ (S7)

where we have abbreviated
$\sum_{\mathbf{r}\mathbf{r}^{\prime},\sigma\sigma^{\prime}}U_{\sigma\sigma^{\prime}}(\mathbf{r},\mathbf{r}^{\prime})\hat{n}_{\mathbf{r}\sigma}\hat{n}_{\mathbf{r}^{\prime}\sigma^{\prime}}\equiv\hat{U}$
and defined the Peierls phase $\Phi_{j}(t)=\mathbf{a}_{j}\cdot\mathbf{A}(t)$.

###  I.2 General Twinning Field

Having obtained the appropriate form for the Hamiltonian, we now consider how
a twinning field might be derived from it. Let us say we have two systems,
described by $\hat{H}^{(1)}$ and $\hat{H}^{(2)}$ each driven by the twinning
vector potential $\mathbf{A}_{O}(t)$, which is designed such that the
expectation of some operator $\hat{O}$ in the two systems satisfies
$O^{(1)}(t)=O^{(2)}(t)$. The Hamiltonian for each system will then be given
by:

$\hat{H}^{(k)}=-t_{0}^{(k)}\sum_{\mathbf{r}^{(k)},\sigma}\sum_{j}\left({\rm
e}^{-i\Phi^{(k)}_{j}(t)}\hat{c}^{\dagger(k)}_{\mathbf{r}^{(k)}\sigma}\hat{c}^{(k)}_{(\mathbf{r}^{(k)}+\mathbf{a}^{(k)}_{j})\sigma}+{\rm
h.c.}\right)+\hat{U}^{(k)},$ (S8)

where $\Phi^{(k)}_{Oj}=\mathbf{a}^{(k)}_{j}\cdot\mathbf{A}_{O}(t)$. In order
to twin $\hat{O}$, it is necessary to match the _derivative_ of that
expectation if the operator itself has no explicit $\mathbf{A}_{O}(t)$
dependence. Given $\frac{{\rm d}O^{(k)}}{{\rm
d}t}=-i\bra{\psi_{k}}\left[\hat{H}^{(k)}(t),\hat{O}(t)\right]\ket{\psi_{k}}$,
for an expectation to be twinned we require

$\bra{\psi_{1}}\left[\hat{H}^{(1)}(t),\hat{O}(t)\right]\ket{\psi_{1}}=\bra{\psi_{2}}\left[\hat{H}^{(2)}(t),\hat{O}(t)\right]\ket{\psi_{2}}.$
(S9)

The conditions under which this equality can be satisfied are best analysed by
breaking up the commutator expectations with the following expectations:

$\displaystyle A_{j}(\psi_{k}){\rm e}^{i\theta_{j}(\psi_{k})}$
$\displaystyle=-t^{(k)}_{0}\sum_{\mathbf{r}^{(k)},\sigma}\braket{\psi_{k}}{\left[\hat{c}^{\dagger(k)}_{\mathbf{r}^{(k)}\sigma}\hat{c}^{(k)}_{(\mathbf{r}^{(k)}+\mathbf{a}^{(k)}_{j})\sigma},\hat{O}\right]}{\psi_{k}},$
(S10) $\displaystyle B^{(k)}(\psi_{k})$
$\displaystyle=\braket{\psi_{k}}{\left[\hat{U}^{(k)},\hat{O}\right]}{\psi_{k}}.$
(S11)

With some rearrangements, we are able to collect all terms dependent on the
driving field, to find a general expression for the twinning field:

$\sum_{j}\left[A_{j}(\psi_{1})\sin\left(\Phi^{(1)}_{Oj}(t)-\theta_{j}(\psi_{1})\right)-A_{j}(\psi_{2})\sin\left(\Phi^{(2)}_{Oj}(t)-\theta_{j}(\psi_{2})\right)\right]=B^{(2)}(\psi_{2})-B^{(1)}(\psi_{1}).$
(S12)

Clearly such an equation only implicitly defines the twinning field and as
such there is no guarantee that a field exists satisfying Eq. (S12). If
however this equation can be reduced to a closed-form expression for the
field, then it is possible to assess when a solution exists.

###  I.3 Twinning Field For Current

While the existence of a twinning field for a general observable depends on if
Eq. (S12) possesses a solution, there is a particular observable for which
twinning fields for a large class of materials can be demonstrated.
Specifically, it is the electronic current $\hat{\mathbf{J}}$, which directly
determines the system’s optical response. This observable is derived directly
from the continuity for the electron density
$\hat{\rho}=\sum_{\mathbf{r},\sigma}\hat{n}_{\mathbf{r}\sigma}$. We obtain an
expression for this operator by evaluating:

$\frac{{\rm d}\hat{n}_{\mathbf{r}\sigma}}{{\rm
d}t}=-i\left[\hat{H},\hat{n}_{\mathbf{r}\sigma}\right]=-it_{0}\left(\sum_{\mathbf{r}^{\prime},j,\sigma^{\prime}}{\rm
e}^{-i\Phi_{j}(t)}\left[\hat{c}^{\dagger}_{\mathbf{r}^{\prime}\sigma^{\prime}}\hat{c}_{(\mathbf{r}^{\prime}+\mathbf{a}_{j})\sigma^{\prime}},\hat{n}_{\mathbf{r}\sigma}\right]-{\rm
h.c.}\right)$ (S13)

Critically, we are able to obtain the second equality due to the fact that the
electron-electron interactions captured by $\hat{U}$ depend only on number
operators (even if higher order terms are included), and such terms will
commute with $\hat{n}_{\mathbf{r}\sigma}$. Evaluating the commutator in the
final equality, we obtain:

$\left[\hat{c}^{\dagger}_{\mathbf{r}^{\prime}\sigma^{\prime}}\hat{c}_{(\mathbf{r}^{\prime}+\mathbf{a}_{j})\sigma^{\prime}},\hat{n}_{\mathbf{r}\sigma}\right]=\delta_{\sigma\sigma^{\prime}}\left(\delta_{\mathbf{r}(\mathbf{r}^{\prime}+\mathbf{a}_{j})}\hat{c}^{\dagger}_{\mathbf{r}^{\prime}\sigma^{\prime}}\hat{c}_{\mathbf{r}\sigma}-\delta_{\mathbf{r}\mathbf{r}^{\prime}}\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{(\mathbf{r}^{\prime}+\mathbf{a}_{j})\sigma^{\prime}}\right).$
(S14)

Defining the operator using ($\left|\mathbf{a}_{j}\right|=a_{j}$)

$\hat{J}_{\mathbf{r}j}=-ia_{j}t_{0}\sum_{\sigma}{\rm
e}^{-i\Phi_{j}(t)}\hat{c}^{\dagger}_{\mathbf{r}\sigma}\hat{c}_{(\mathbf{r}^{\prime}+\mathbf{a}_{j})\sigma}-{\rm
h.c.}$ (S15)

using this operator, we find that the continuity equation can be expressed as

$\frac{{\rm d}\hat{\rho}}{{\rm
d}t}=\sum_{j}\frac{1}{a_{j}}\sum_{\mathbf{r}}\left(\hat{J}_{\mathbf{r}j}-\hat{J}_{(\mathbf{r}-\mathbf{a}_{j})j}\right).$
(S16)

The right hand side of this equation can immediately be recognised as a
discretised divergence, and hence we can identify the current in the $j$
direction as $\hat{J}_{j}=\sum_{\mathbf{r}}\hat{J}_{\mathbf{r}j}$.

In the case of this observable the expression for the twinning field
($\Phi^{(k)}_{\tau j}(t)=\mathbf{a}^{(k)}_{j}\cdot\mathbf{A}_{\tau}(t)$) is
significantly simplified. By first expressing the nearest neighbour hopping
expectation in a polar form

$\sum_{\mathbf{r}^{(k)},\sigma}\braket{\psi_{k}}{\hat{c}^{\dagger(k)}_{\mathbf{r}^{(k)}\sigma}\hat{c}^{(k)}_{(\mathbf{r}^{(k)}+\mathbf{a}^{(k)}_{j})\sigma}}{\psi_{k}}=R_{j}(\psi_{k}){\rm
e}^{i\theta_{j}(\psi_{k})}$ (S17)

The current expectation in the $j$ direction may be given as

$J^{(k)}_{j}(t)=\braket{\psi_{k}}{\hat{J}_{j}}{\psi_{k}}=2a^{(k)}_{j}t^{(k)}_{0}R_{j}(\psi_{k})\sin(\Phi^{(k)}_{\tau
j}(t)-\theta_{j}(\psi_{k})).$ (S18)

Equating $\hat{J}_{j}^{(1)}=\hat{J}_{j}^{(2)}$, we obtain

$a^{(1)}_{j}t^{(1)}_{0}R_{j}(\psi_{1})\sin(\Phi^{(1)}_{\tau
j}(t)-\theta_{j}(\psi_{1}))-a^{(2)}_{j}t^{(2)}_{0}R_{j}(\psi_{2})\sin(\Phi^{(2)}_{\tau
j}(t)-\theta_{j}(\psi_{2}))=0.$ (S19)

Comparing to Eq. (S12), we find two major simplifications that bring us close
to a closed form equation for the field. First, the twinning field does not
explicitly depend on any expectation of the potential $\hat{U}^{(k)}$ term
(although the evolution of the state will depend on this potential). Secondly,
the specific form of the current operator means that each component of the
current in the two systems can be equated separately, meaning one has an
equation for the twinning field in each of the $j$ nearest neighbour
directions, rather than a single equation summing over all components.

Unfortunately, it is still not possible to obtain a closed form equation for
the twinning field from Eq. (S19). To amend this we make the further
assumption that both systems have the same lattice structure,
$\mathbf{r}^{(1)}=\mathbf{r}^{(2)}$ and
$\mathbf{a}^{(1)}_{j}=\mathbf{a}^{(2)}_{j}$. In this case $\Phi^{(1)}_{\tau
j}(t)=\Phi^{(2)}_{\tau j}(t)\equiv\Phi_{\tau j}(t)$, and Eq. (S19) can be
rearranged to express the twinning field directly in terms of the systems’
expectations:

$\displaystyle\Phi_{\tau j}(t)$
$\displaystyle=\arctan\left(\xi_{j}(\psi_{1},\psi_{2})\right)$ (S20)
$\displaystyle\xi_{j}(\psi_{1},\psi_{2})$ $\displaystyle=\frac{\lambda
R_{j}(\psi_{1})\sin(\theta_{j}(\psi_{1}))-R_{j}(\psi_{2})\sin(\theta_{j}(\psi_{2}))}{\lambda
R_{j}(\psi_{1})\cos(\theta_{j}(\psi_{1}))-R_{j}(\psi_{2})\cos(\theta_{j}(\psi_{2}))},$
(S21)

using $\lambda=\frac{t^{(1)}_{0}}{t^{(2)}_{0}}$. Note that both
$R_{j}(\psi_{k})$ and $\theta_{j}(\psi_{k})$ will have definite values for
every valid state $\ket{\psi_{k}}$ in each system. If this were not the case
then the action of the nearest neighbour hopping operator (and therefore the
Hamiltonian) on the system would not be well defined. For this reason,
$\xi_{j}(\psi_{1},\psi_{2})$ will have a definite real value regardless of the
states of the system. Given the $\arctan$ function maps any real to
$(-\pi/2,+\pi/2)$, we conclude that $\Phi_{\tau j}$ will always have a
definite value.

To summarise, while an expression for the twinning field can always be
obtained, in general it cannot be ascertained as to whether this equation will
possess a solution. _If the two systems share an identical lattice structure
then it it possible to express the twinning field in a closed form._ Under
these conditions, the twinning field may be explicitly expressed as a function
of well-defined system expectations, and given that the twinning field is a
function whose range and domain extends over the entire real line, it is
guaranteed to exist regardless of the value of those system expectations. This
leads to the perhaps surprising conclusion that ultimately it is the lattice
structure of a material, rather than its electronic interactions, that
determine if its optical response can be twinned.

###  I.4 Specification to Hexagonal Bipartite Lattice

Here we outline the notation used for the twinning field for a hexagonal
bipartite lattice (shown in figure 3), in order to simulate graphene and hBN.
In this case, the bipartite nature of the lattice means that the sum over
lattice positions can be indexed directly, where sites belonging to the $A$
and $B$ sublattices are contained in the sets $\mathcal{A}$ and $\mathcal{B}$
respectively. If an onsite potential is used, the Hamiltonian reads:

$\hat{H}=-t_{0}\sum_{<i,j>,\sigma}\left({\rm
e}^{-i\Phi_{j}(t)}\hat{c}_{i\sigma}^{\dagger}\hat{b}_{j\sigma}+{\rm
h.c.}\right)+\sum_{i\in\mathcal{A}}U_{A}\hat{n}_{i\uparrow}^{A}\hat{n}_{i\downarrow}^{A}+\sum_{j\in\mathcal{B}}U_{B}\hat{n}_{j\uparrow}^{B}\hat{n}_{j\downarrow}^{B}$
(S22)

here $\hat{c}$ operators act on the $A$ sublattice, while $\hat{b}$ acts on
the $B$ sublattice. Each sublattice possesses a separate onsite potential
$U_{A}$ and $U_{B}$, using $\hat{n}_{i\sigma}^{A/B}$ as the number operators
acting on the appropriate lattice. Finally, each site has three neighbours
related by the vectors:

$\displaystyle\mathbf{a}_{1}=a\begin{bmatrix}1\\\ 0\end{bmatrix},\ \
\mathbf{a}_{2}=\frac{a}{2}\begin{bmatrix}1\\\ -\sqrt{3}\end{bmatrix},\ \
\mathbf{a}_{3}=\frac{a}{2}\begin{bmatrix}1\\\ \sqrt{3}\end{bmatrix}.$ (S23)

A current in each of these directions can be defined in the same manner as the
previous section:

$\hat{J}_{j}=-iat_{0}\sum_{i\in\mathcal{A},\sigma}\left({\rm
e}^{-i\Phi_{j}(t)}\hat{c}_{i\sigma}^{\dagger}\hat{b}_{j\sigma}-{\rm
h.c.}\right).$ (S24)

Then, if we once again define the functionals for each separate system,

$\sum_{i\in\mathcal{A}^{(k)},\sigma}\braket{\psi_{k}}{\hat{c}^{\dagger(k)}_{i\sigma}\hat{b}^{(k)}_{j\sigma}}{\psi_{k}}=R_{j}(\psi_{k}){\rm
e}^{i\theta_{j}(\psi_{k})},$ (S25)

it is possible to apply Eq. (S20) directly.

Even in this case (where each component of a twinning field has a closed
form), it is important to remark that in general this does not mean the fields
generated are physically realisable. To see this, consider the fact that one
has three independent equations for the $\Phi_{j}(t)$, but geometrically there
are only two independent components. In Cartesian coordinates, the fields in
the three lattice directions are given by:

$\Phi_{\tau 1}(t)=\Phi_{\tau x}(t),\ \ \Phi_{\tau
2}(t)=\frac{1}{2}\left(\Phi_{\tau x}(t)+\sqrt{3}\Phi_{\tau y}(t)\right),\ \
\Phi_{\tau 3}(t)=\frac{1}{2}\left(\Phi_{\tau x}(t)-\sqrt{3}\Phi_{\tau
y}\right).$ (S26)

While it is always possible to find a solution for the fields in the lattice
directions, there is no guarantee that these will be consistent with the
geometric composition of these fields. Fortunately, when the field is
polarised along one of the lattice directions (say $\mathbf{a}_{1}$), the
symmetry of the system guarantees that currents in the two remaining
directions will be equivalent, i.e. $J_{2}(t)=J_{3}(t)$, such that
$J_{y}(t)=0$. Furthermore, local current conservation means that in this case
$J_{1}(t)=J_{2}(t)+J_{3}(t)$. Consequently, the total current in the $x$
direction will be $J_{x}(t)=2J_{1}(t)$. This means only the component of the
twinning field in that direction need be calculated, removing the
aforementioned geometric constraints.

###  I.5 Uniqueness of Twinning Fields

While a solution for the twinning field may exist, the existence of this
solution says nothing about the _number_ of twinning fields that could exist
for a given pair of systems – i.e. is a twinning field unique? It is possible
to deduce an answer to this question by substituting $\Phi_{\tau}(t)$ into the
Hamiltonian for each system. In this case we obtain an explicitly nonlinear
evolution governed by a ‘twinning Hamiltonian’ for the joint system (here
using a 1D example):

$\displaystyle\hat{H}^{\rm tot}_{\tau}(\psi_{1},\psi_{2})$
$\displaystyle=\hat{H}_{\tau}^{(1)}(\psi_{1},\psi_{2})\otimes\mathbb{1}^{(2)}+\mathbb{1}^{(1)}\otimes\hat{H}_{\tau}^{(2)}(\psi_{1},\psi_{2})$
(S27) $\displaystyle\hat{H}^{(k)}_{\tau}\left(\psi_{1},\psi_{2}\right)$
$\displaystyle=\sum_{\sigma,j}\left[P^{(k)}_{+}\hat{c}_{j\sigma}^{\dagger}\hat{c}_{j+1\sigma}+P^{(k)}_{-}\hat{c}_{j\sigma}\hat{c}^{\dagger}_{j+1\sigma}\right]+\hat{U}^{(k)},$
(S28) $\displaystyle P^{(k)}_{\pm}$
$\displaystyle=-\frac{t^{(k)}_{0}}{\sqrt{1+\xi^{2}(\psi_{1},\psi_{2})}}\left(1\pm
i\xi(\psi_{1},\psi_{2})\right).$ (S29)

Using this twinning Hamiltonian, we can assess the conditions under which the
twinning field is unique. If the evolution generated by the Hamiltonian in Eq.
(S27) from a given initial condition has only one solution, it implies there
is only one field which satisfies the twinning condition. Applying Lipschitz
analysis 111In order to assess the uniqueness of the evolution of the joint
state $\ket{\Psi}=\ket{\psi_{1}}\otimes\ket{\psi_{2}}$ under $\hat{H}^{\rm
tot}_{\tau}$, we make use of the Picard-Lindelöf theorem [4]. This states that
the solution to a differential equation
$y^{\prime}(t)=f(t,y(t)),y(t_{0})=y_{0}$ is unique if $f(t,y(t))$ is a
_Lipschitz continuous_ (LC) [5] function over all $y$. It therefore follows
that the conditions under which $\Phi_{\tau}(t)$ is unique are equivalent to
those for which $\hat{H}^{\rm tot}_{\tau}\ket{\Psi}$ is LC. to this
Hamiltonian, it is possible to prove the following:

_For a finite dimensional system (i.e., when the twinning Hamiltonian
$\hat{H}^{\rm tot}_{\tau}$ is a finite dimensional matrix), if $\ket{\Psi}$
solves Eq. (S27) and satisfies the constraints_

$\displaystyle\left|R(\psi_{1/2})\right|>$ $\displaystyle\epsilon_{1}>0,$
(S30) $\displaystyle\left|{\rm Re}\left(\lambda
K(\psi_{1})-K(\psi_{2})\right)\right|>$ $\displaystyle\epsilon_{2}>0,$ (S31)

_where $\epsilon_{1/2}$ are any constants, then the twinning field
$\Phi_{\tau}(t)$ is unique._

The first condition [Eq. (S30)] can be straightforwardly interpreted as the
requirement that there be some kinetic energy in each system, a natural
precondition for observing a current. Some sense can be made of the second
condition by noting that evolution under the twinning field implies
$\frac{{\rm d}J^{(1)}(t)}{{\rm d}t}=\frac{{\rm d}J^{(2)}(t)}{{\rm d}t}$.
Evaluating this derivative, we obtain:

$\displaystyle\dot{\Phi}_{\tau}(t){\rm Re}\left(\lambda
K(\psi_{1})-K(\psi_{2})\right)={\rm Re}$
$\displaystyle\left(\left\langle\psi_{1}(t)\left|\left[\hat{U}^{(1)},\hat{K}\right]\right|\psi_{1}(t)\right\rangle\right)-{\rm
Re}\left(\left\langle\psi_{2}(t)\left|\left[\hat{U}^{(2)},\hat{K}\right]\right|\psi_{2}(t)\right\rangle\right).$
(S32)

Clearly, when Eq. (S31) is violated, $\dot{\Phi}_{\tau}$ can take any value
while still functioning as a twinning field. This means that in principle, if
a twinning field is non-unique, there are an _infinite_ number of solutions.
These solutions will however be constrained to follow the same trajectory up
to the point where one of Eq. (S30) or Eq. (S31) is violated, and a
proliferation of possible trajectories occurs. This non-uniqueness is
therefore qualitatively different to that previously observed in control
fields for tracking control [2, 3], where a branch point generates a well-
defined number of additional solutions.

## References

  * Note [1] In order to assess the uniqueness of the evolution of the joint state $\mathinner{|{\Psi}\rangle}=\mathinner{|{\psi_{1}}\rangle}\otimes\mathinner{|{\psi_{2}}\rangle}$ under $\hat{H}^{\rm tot}_{\tau}$, we make use of the Picard-Lindelöf theorem [4]. This states that the solution to a differential equation $y^{\prime}(t)=f(t,y(t)),y(t_{0})=y_{0}$ is unique if $f(t,y(t))$ is a _Lipschitz continuous_ (LC) [5] function over all $y$. It therefore follows that the conditions under which $\Phi_{\tau}(t)$ is unique are equivalent to those for which $\hat{H}^{\rm tot}_{\tau}\mathinner{|{\Psi}\rangle}$ is LC.
  * McCaul  _et al._ [2020] G. McCaul, C. Orthodoxou, K. Jacobs, G. H. Booth, and D. I. Bondar, Controlling arbitrary observables in correlated many-body systems, Phys. Rev. A 101, 053408 (2020).
  * Jha  _et al._ [2009] A. Jha, V. Beltrani, C. Rosenthal, and H. Rabitz, Multiple solutions in the tracking control of quantum systems, J. Phys. Chem. A 113, 7667 (2009).
  * Nagle [2011] R. K. Nagle,  _Fundamentals of Differential Equations and Boundary Value Problems (6th Edition) (Featured Titles for Differential Equations)_ (Pearson, 2011).
  * Folland [2007] G. B. Folland,  _Real Analysis: Modern Techniques and Their Applications_ (Wiley, 2007).


[b]J. T. Kuhlmann

# On improvement of the axial-vector current with stabilised Wilson fermions

P. Fritzsch  J. Heitger

###### Abstract

We report on the determination of the improvement coefficient $c_{\rm A}$ for
the non-singlet axial-vector current $A^{a}_{\mu}(x)$ in the framework of
stabilised Wilson-Clover fermions within three flavour lattice QCD. This is
done by requiring the PCAC relation to hold for two different pseudo-scalar
states. To generate these states, wavefunctions altering spatial structures on
the boundaries of a Schrödinger functional lattice are employed and some
variations of the (previously applied) wavefunction method are explored. The
improvement coefficient is determined on a few ensembles for a range of gauge
couplings that is potentially useful for future applications. Preliminary
results on the renormalisation constants $Z_{\rm V}$ for the vector current
and $Z_{\rm A}$ for the axial-vector current are also presented.

##  1 Introduction

Leptonic pseudo-scalar hadronic decays have a significant impact on precisely
estimated entries of the CKM matrix and therefore are of great interest for
indirect tests of the Standard Model. Since these decays involve the QCD
hadronic matrix element of the axial-vector current

$\langle
0|A^{a}_{\mu}(x)|\pi^{b}(p)\rangle=-ip_{\mu}f_{\pi}\,\delta^{ab}e^{-ipx}\,,$
(1)

lattice QCD is the favoured way to reliably estimate the low-energy part of
the theory contribution to these tests.

As is well-known, Wilson fermions, while having many advantages, break chiral
symmetry and need to be improved to remove its leading ${\rm O}(a)$
discretisation errors. For this purpose, resorting to the well-established
Symanzik improvement programme [1], the fermion action receives an extra term
known as the Sheikholeslami-Wohlert term. Quark bilinears also need to be
improved by adding suitable extra terms. Within this programme it turns out
that the improvement term of the axial-vector current is proportional to the
derivative of the pseudo-scalar density. In a mass independent scheme, the
(multiplicatively) renormalised and improved current for degenerate quarks is
of the form

$(A_{{\rm R}})^{a}_{\mu}(x)=Z_{\rm A}(1+b_{\rm A}m_{q})(A_{{\rm
I}})^{a}_{\mu}(x)=Z_{\rm A}(1+b_{\rm A}m_{q})\left[A^{a}_{\mu}(x)+c_{\rm
A}\partial_{\mu}P^{a}(x)\right]$ (2)

with the axial-vector current and pseudo-scalar density being defined as

$A^{a}_{\mu}\left(x\right)=\bar{\psi}\left(x\right)\gamma_{\mu}\gamma_{5}\frac{\lambda^{a}}{2}\psi\left(x\right)\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
}P^{a}=\bar{\psi}\left(x\right)\gamma_{5}\frac{\lambda^{a}}{2}\psi\left(x\right)\,.$
(3)

In practice, for a non-perturbative determination of $c_{\rm A}$ and $Z_{\rm
A}$ one works in the chiral limit of vanishing quark masses where the $b_{\rm
A}$-term vanishes.

Here we report on our efforts to improve and renormalise the axial-vector
current for three-flavour lattice QCD with stabilised Wilson fermions (SWF)
[2]. For this computation we employ ensembles with fairly large volumes
($L\approx 3\,\rm{fm}$), Schrödinger functional boundary conditions and
lattice constants in a range between $a=0.12\,\rm{fm}$ ($\beta=3.685$) and
$a=0.055\,\rm{fm}$ ($\beta=4.10$) in the chiral limit as well as at the
symmetric point of fully degenerate sea and valence quark masses. Since SWF
provide a new ${\rm O}(a)$ discretisation of QCD, their improvement
coefficients and renormalisation constants differ from the known $N_{\rm f}=3$
values (see e.g. [3, 4, 5, 6]) for ordinary Wilson-Clover fermions with the
same, tree-level Symanzik-improved gauge action.

##  2 Stabilised Wilson formulation

Numerical simulations with the Wilson-Clover Dirac operator

$D_{{\rm
W}}=\frac{1}{2}\left\\{\gamma_{\mu}\left(\nabla_{\mu}^{*}+\nabla_{\mu}\right)-\nabla_{\mu}^{*}\nabla_{\mu}\right\\}+c_{\rm
sw}\frac{i}{4}\sigma_{\mu\nu}F_{\mu\nu}+m_{0}$ (4)

require stabilising measures as proposed in [2]. Based on the observation that
the clover term in $D_{{\rm W}}$ hinders stable inversion of the Dirac matrix,
an alternative form of the Dirac operator was introduced, in which the on-site
terms are exponentiated, i.e.

$M_{0}+c_{\rm sw}\frac{i}{4}\sigma_{\mu\nu}F_{\mu\nu}\quad\rightarrow\quad
M_{0}\exp\left(\frac{c_{\rm
sw}}{M_{0}}\frac{i}{4}\sigma_{\mu\nu}F_{\mu\nu}\right)\,,$ (5)

with $M_{0}=m_{0}+4$. The exponential ensures an intrinsic bound from below in
this term. As a consequence, the Dirac operator does receive fewer near-zero
eigenvalues, stabilising the calculation. This approach was originally
developed for master-field simulations [2], but it turns out to also be
beneficial for light quarks, coarse gauge fields and large lattices [7]. This
is the situation we are focussing on.

While the Sheikholeslami-Wohlert coefficient $c_{\rm sw}$ for the new action
has already been determined for a range of couplings [2, 7], the improvement
coefficients for quark-bilinears such as the axial-vector current are yet to
be computed.

##  3 Improvement of the axial-vector current

Along the lines of [8, 4], using the Schrödinger functional framework, one can
derive an expression for the coefficient $c_{\rm A}$. To do so, we introduce
the correlation functions

$f_{\rm A}(x_{0};\omega)=-\frac{a^{3}}{3L^{6}}\sum_{\vec{x}}\langle
A^{a}_{0}(x)O^{a}(\omega)\rangle\text{\quad and\quad}f_{\rm
P}(x_{0};\omega)=-\frac{a^{3}}{3L^{6}}\sum_{\vec{x}}\langle
P^{a}(x)O^{a}(\omega)\rangle$ (6)

with the boundary source operator at $x_{0}=0$ being defined as

$O^{e}(\omega)=a^{6}\sum_{\vec{y},\vec{z}}\bar{\zeta}(\vec{y})\gamma_{5}\frac{\lambda^{e}}{2}\omega(\vec{y}-\vec{z})\zeta(\vec{z})\,,$
(7)

where $\zeta$ is the boundary quark field. Similar correlation functions can
also be introduced with sources at $x_{0}=T$. These are built from boundary
operators $O^{\prime f}(\omega)$ and complementary boundary quark fields
$\zeta^{\prime}$. In this context, $\omega$ denotes an arbitrary spatial
structure on the boundary that we will call wavefunction. The first step in
our analysis is to find suitable wavefunctions such that the ground and first
excited state can be readily extracted from the correlators. These are
approximated as in [8, 4] by constructing linear combinations of basis
wavefunctions. Here we consider a set of five functions of the form

$\omega^{\text{b}1}=\mathrm{e}^{-\frac{r}{a_{0}}}\,,\leavevmode\nobreak\
\omega^{\text{b}2}=r\leavevmode\nobreak\
\mathrm{e}^{-\frac{r}{a_{0}}}\,,\leavevmode\nobreak\
\omega^{\text{b}3}=\mathrm{e}^{-\frac{r}{2a_{0}}}\,,\leavevmode\nobreak\
\omega^{\text{b}4}=\text{const}\,,\leavevmode\nobreak\
\omega^{\text{b}5}=r^{2}\leavevmode\nobreak\ \mathrm{e}^{-\frac{r}{a_{0}}}\,,$
(8)

which resemble those of the hydrogen atom, with $r=|\vec{x}-\vec{y}|$ and
$a_{0}$ parametrising the spatial extent of the function. To find a well-
suited linear combination, the boundary-to-boundary correlator,

$F_{1}(\omega^{{\rm b}i},\omega^{{\rm b}\\!j})=-\frac{1}{3L^{6}}\langle
O^{\prime a}(\omega^{{\rm b}\\!j})O^{a}(\omega^{{\rm b}i})\rangle\,,$ (9)

which for our purposes can be viewed as a $5\\!\times\\!5$-matrix in
wavefunction space, is evaluated. We can pick any subset of indices as a
submatrix111These need to be at least $3\\!\times\\!3$-submatrices for the
eigenvectors to be fully independent of each other. and calculate its
eigenvectors, sorted by eigenvalue. The „standard“ choice that was previously
used in refs. [8, 4] can be regained by setting $i,j\in\\{1,2,3\\}$. In our
setup it turned out that with this choice $c_{\rm A}$ inherits some ambiguity
that may be avoided by other linear combinations. Hence we also explored the
submatrix with $i,j\in\\{1,2,4\\}$. Exploiting the overlap information of the
basis wavefunctions to approximately prepare the ground and first excited
state then amounts to projecting the correlators $f_{\rm A}$ and $f_{\rm P}$
with given basis wavefunctions onto the eigenvectors belonging to the largest
two eigenvalues.

The improvement coefficient $c_{\rm A}$ is determined in the following way:
The starting point is the PCAC relation

$\partial_{\mu}\langle A^{a}_{\mu}(x)O^{a}\rangle=2m_{\rm
PCAC}\leavevmode\nobreak\ \langle P^{a}(x)O^{a}\rangle\,,$ (10)

which follows from an invariance under chiral rotation and holds as an
operator identity in the continuum theory. Therefore, all effects that break
it on the lattice must stem from the discretisation. Requiring the identity to
be satisfied up to ${\rm O}(a^{2})$ cut-off effects amounts to inserting the
improved current from eq. (2) into this relation and thus yields for the ${\rm
O}(a)$ improved PCAC quark mass

$am_{\rm PCAC}\equiv\frac{a\leavevmode\nobreak\ \partial_{\mu}\langle(A_{{\rm
I}})^{a}_{\mu}(x)O^{a}\rangle}{2\leavevmode\nobreak\ \langle
P^{a}(x)O^{a}\rangle}+{\rm O}(a^{2})\,,$ (11)

provided $c_{\rm A}$ is fixed non-perturbatively. As detailed in the previous
works, it can be rewritten in terms of the Schrödinger functional correlation
functions from above as

$a\leavevmode\nobreak\ m_{\rm
PCAC}(x_{0};\omega^{i})=\frac{a\leavevmode\nobreak\ \hat{\partial}_{0}f_{\rm
A}(x_{0};\omega^{i})+c_{\rm A}\leavevmode\nobreak\ a^{2}\leavevmode\nobreak\
\hat{\partial}^{2}_{0}f_{\rm P}(x_{0};\omega^{i})}{2\leavevmode\nobreak\
f_{\rm P}(x_{0};\omega^{i})}\,.$ (12)

Decomposing the PCAC mass as

$m_{\rm PCAC}(x_{0};\omega^{i})=r(x_{0};\omega^{i})+c_{\rm
A}\leavevmode\nobreak\ a\leavevmode\nobreak\ s(x_{0};\omega^{i})\,,$ (13)

with

$r(x_{0};\omega^{i})=\frac{\hat{\partial}_{0}f_{\rm
A}(x_{0};\omega^{i})}{2\leavevmode\nobreak\ f_{\rm
P}(x_{0};\omega^{i})}\text{\quad
and\quad}s(x_{0};\omega^{i})=\frac{\hat{\partial}^{2}_{0}f_{\rm
P}(x_{0};\omega^{i})}{2\leavevmode\nobreak\ f_{\rm P}(x_{0};\omega^{i})}\,,$
(14)

and assuming that it holds for the two lowest states in the pseudo-scalar
channel (isolated as described above) labelled by indices $i=0,1$, we can
solve for $c_{\rm A}$:

$c_{\rm
A}(x_{0})=-\frac{r(x_{0};\omega^{1})-r(x_{0};\omega^{0})}{a\leavevmode\nobreak\
(s(x_{0};\omega^{1})-s(x_{0};\omega^{0}))}=-\frac{\Delta
r(x_{0})}{a\leavevmode\nobreak\ \Delta s(x_{0})}\,.$ (15)

Note that, as observed in [8, 4], the precision of the estimate for $c_{\rm
A}$ is dominated by the statistical error of $a\Delta s$. Therefore an
examination of the $x_{0}$-dependence of the functions $\Delta r$ and $a\Delta
s$ may be helpful. As one expects $c_{\rm A}(x_{0})$ in eq. (15) to develop a
plateau over a certain range of $x_{0}$, a suitable plateau range needs to be
identified, in which the noise in $\Delta r$ and $a\Delta s$ is small.

##  4 Renormalisation

For the calculation of the renormalisation factor $Z_{\rm A}$ we adopt the
strategy explained in [9, 5]. We only recall here that the calculation
involves a Ward identity originating from an infinitesimal chiral rotation of
the fields entering the expectation value

$\langle
A_{\nu}^{b}(y)O_{\text{ext}}\rangle=\int\mathcal{D}[\bar{\psi},\psi,U]A^{b}_{\nu}(y)O_{\text{ext}}e^{-S}\,,$
(16)

where we include the current into the chiral rotation and obtain from the
variational principle

$0=\langle\delta A^{b}_{\nu}(y)O_{\text{ext}}\rangle-\langle\delta
SA^{b}_{\nu}(y)O_{\text{ext}}\rangle\,.$ (17)

For $O_{\text{ext}}$ which we set

$O_{\text{ext}}=O^{ef}_{\text{ext}}=O^{e}(O^{\prime})^{f}\,,$ (18)

where $O^{e}$ and $(O^{\prime})^{f}$ are the pseudo-scalar boundary sources as
defined in eq. (7).

Omitting all intermediate steps, translating the expressions into their
lattice counterparts and formulating them in terms of renormalised, improved
Schrödinger functional correlators, one can finally isolate $Z_{\rm A}$ as

$Z_{\rm A}=\sqrt{\frac{F_{1}}{F_{\rm{AA}}^{\rm I}(x_{0},y_{0})-2m_{\rm
PCAC}\tilde{F}_{\rm{PA}}^{\rm I}(x_{0},y_{0})}}\,,$ (19)

where $F_{\rm{AA}}^{\rm I}$ is defined as

$F_{\rm{AA}}^{\rm I}=F_{\rm{AA}}(x_{0},y_{0})+ac_{\rm
A}\left(\partial_{x_{0}}F_{\rm{PA}}(x_{0},y_{0})+\partial_{y_{0}}F_{\rm{AP}}(x_{0},y_{0})\right)+a^{2}c_{\rm
A}^{2}\left(\partial_{x_{0}}\partial_{y_{0}}F_{\rm{PP}}(x_{0},y_{0})\right)$
(20)

and $\tilde{F}_{\rm{PA}}^{\rm{I}}$ as

$\tilde{F}^{\rm
I}_{\rm{PA}}(x_{0},y_{0})=\tilde{F}_{\rm{PA}}(x_{0},y_{0})+c_{\rm
A}\partial_{y_{0}}\tilde{F}_{\rm{PP}}(x_{0},y_{0})\,,$ (21)

with

$F_{\rm{XY}}(x_{0},y_{0})\equiv-2a^{6}\sum_{\vec{x},\vec{y}}\langle
X^{1}(x)Y^{2}(y)O^{2}(O^{\prime})^{1}\rangle$ (22)

and

$\tilde{F}_{\rm{XY}}(x_{0},y_{0})\equiv
a\sum_{z_{0}=x_{0}}^{y_{0}}w(z_{0})F_{\rm{XY}}(z_{0},y_{0})\text{\quad
with\quad}w(z_{0})=\begin{cases}\frac{1}{2}\;\text{if\leavevmode\nobreak\
}z_{0}=x_{0}\text{\;or\;}z_{0}=y_{0}\\\ 1\;\text{else}\end{cases}\,.$ (23)

Each of the 4-point functions receives contributions from six different
diagrams shown in fig. 1.

Figure 1: The six non-vanishing Wick contractions contributing to the 4-point
functions $F_{{\rm XY}}$.

Two of these are disconnected and can be evaluated as a product of two
separate 2-point functions.

Additionally, in this calculation we consider the vector current which enters
the following Ward identity:

$\int\mathrm{d}^{3}y\langle
V^{c}_{0}(y)O^{e}(O^{\prime})^{f}\rangle=if^{ce}_{g}\langle
O^{g}(O^{\prime})^{f}\rangle\,,$ (24)

where $O^{e}$ and $(O^{\prime})^{f}$ are defined as before. From this, a
simple expression for $Z_{\rm V}$ can be derived [10]:

$Z_{\rm V}=\frac{F_{1}}{f_{\rm V}}+{\rm O}(a^{2})\,,$ (25)

in which $f_{\rm V}$ is

$f_{\rm V}(x_{0})=i\sum_{x}\varepsilon_{abc}\langle O^{\prime
a}V^{b}_{0}(x_{0})O^{c}\rangle\,.$ (26)

##  5 Ensembles and error analysis

After having discussed the theoretical aspects, we introduce the ensembles
used in the numerical computations so far. To make our results on $c_{\rm A}$,
$Z_{\rm A}$ and $Z_{\rm V}$ useful for other groups, the range of ensembles
covered are at the same lattice spacings $a$ as those generated by
OpenLat222see also: https://openlat1.gitlab.io/ [7], only differing in
boundary conditions. The ensembles used here are listed in tab. 1. Our error
analysis is done employing the $\Gamma$-method [11] using the python-
implementation pyerrors described in [12].

Table 1: Ensembles with Schrödinger functional boundary conditions used in this project. The mass point indicates to which line of constant physics the ensemble belongs. $T$ |  | $L^{3}$ |  | $\beta$ |  | $\kappa_{u/d/s}$ |  |  $a$ [fm] |  | mass point |  | # of configs |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$24$ |  | $16^{3}$ |  | $3.80$ |  | $0.1392500$ |  | $0.095$ |  | chir. |  | $2875$ |
$24$ |  | $24^{3}$ |  | $3.80$ |  | $0.1392500$ |  | $0.095$ |  | chir. |  | $2523$ |
$32$ |  | $32^{3}$ |  | $3.80$ |  | $0.1392500$ |  | $0.095$ |  | chir. |  | $1167$ |
$24$ |  | $24^{3}$ |  | $3.685$ |  | $0.1394400$ |  | $0.120$ |  | symm. |  | $540$ |
$32$ |  | $32^{3}$ |  | $3.80$ |  | $0.1389630$ |  | $0.095$ |  | symm. |  | $956$ |
$56$ |  | $56^{3}$ |  | $4.10$ |  | $0.1380000$ |  | $0.055$ |  | symm. |  | $33$ |

As stated in the introduction, we are interested in two lines of constant
physics: at the chiral point of (almost) vanishing quark masses and at the
symmetric point of fully degenerate massive quark flavours.

##  6 Results

|
---|---
Figure 2: Preliminary results for $c_{\rm A}$ at $\beta=3.80$ ($a\approx
0.094\,\text{fm}$, left) at the chiral point and $\beta=4.10$ ($a\approx
0.055\,\text{fm}$, right) at the symmetric point with the standard projection.
Figure 3: Example for the possible projections on the $32^{4}$ lattice at
$\beta=3.80$ at the chiral point. Data points shown here are the plateau-
values of $x_{0}\in[T/4-2,T/4+2]$. The boxes indicate the spread of $c_{\rm
A}$ in the plateau region. The standard projection is marked in red, while the
new, preferred projection is marked in green.

In the analysis for $c_{\rm A}$ we have seen that the quality and longevity of
the plateau may be problematic in the large-volume setup adopted here. For the
preliminary results we set the plateau region as $x_{0}\in[T/4-2,T/4+2]$,
which works out fine in most cases, as can be seen for the two representative
examples in fig. 2.

As mentioned above, it has also been seen in our tests that the standard
projection (used in [4, 8]) is not always optimal. An alternative was found by
scanning all possible projections which can be built after diagonalising
$F_{1}$ on various (subsets of) basis wavefunctions. The outcome of these
scans can be inferred from fig. 3. While the standard projection shown in red,
assessed in isolation, yields a legitimate value, the majority of the other
projections are placed systematically around a different value. This could
hint at an overestimation of $c_{\rm A}$ (in magnitude) when using the
standard projection. However, the new, preferred projection exhibits smaller
errors and yields results consistent with most other projections, thus
appearing to be less affected by systematics due to the choice of the
wavefunction basis. The spread of $c_{\rm A}$ in the plateau region (shown in
fig. 3 by the vertical extent of the boxes) is also smaller for this
projection, which endorses the better quality of the associated plateau.

In fig. 4 a preliminary interpolation in $g_{0}^{2}$ of the results on the
symmetric point ensembles is displayed. So far, the ensembles studied do not
yet allow for an analogous interpolation in the chiral case. However, we can
see that in comparison with the results for ordinary Wilson-Clover fermions
[4] the estimates for $c_{\rm A}$ appear to be smaller in magnitude which may
hint at smaller cut-off effects for stabilised Wilson fermions at the same
lattice spacing. Of course, whether this is a general feature of this
discretisation still needs to be confirmed in phenomenological applications.

For the renormalisation constants, preliminary results are listed in tables. 2
and 3. These indicate that the estimates obtained are within the same ballpark
as the ones for ordinary Wilson-Clover fermions in $N_{\rm f}=3$ [5, 13].

|
---|---
Figure 4: Preliminary $g_{0}^{2}$-dependence of $c_{\rm A}$ for ensembles at the symmetric point with $\beta=3.685$, $\beta=3.80$, $\beta=4.10$ with the standard (left) and the preferred (right) projection. In orange, the point for the chiral case at $\beta=3.80$ is shown. The red point gives $N_{\rm f}=3$ result for traditional Wilson-Clover fermions from [4]. Table 2: Preliminary estimates of renormalisation constants for the axial-vector channel. The labels f, c and m refer to full, connected and massive definitions of $Z_{\rm A}$; see [5] for details. Ensemble | $Z_{\rm A}^{{\rm f}}$ | $Z_{\rm A}^{\rm f,m}$ | $Z_{\rm A}^{{\rm c}}$ | $Z_{\rm A}^{\rm c,m}$ |
---|---|---|---|---|---
$24\\!\times\\!16^{3}$, $\beta=3.80$, chir. | $0.7638(76)$ | $0.7585(80)$ | $0.753(20)$ | $0.773(21)$ |
$24\\!\times\\!24^{3}$, $\beta=3.80$, chir. | $0.779(15)$ | $0.741(13)$ | $0.7631(50)$ | $0.7782(39)$ |

As the 4-point functions used here are prone to large statistical
fluctuations, it seems more promising to determine $Z_{\rm A}$ in the chirally
rotated Schrödinger functional ($\chi\rm{SF}$) scheme [6, 14], because in this
framework only 2-point functions are involved.

Figure 5: Example for the calculation of $Z_{\rm V}$ on the chiral ensemble at
$\beta=3.80$ ($a\approx 0.094\,\text{fm}$).

Table 3: Preliminary estimates of renormalisation constants for the vector
channel.

Ensemble | $Z_{\rm V}$
---|---
$24\\!\times\\!16^{3}$, $\beta=3.80$, chir. | $0.729141(94)$
$24\\!\times\\!24^{3}$, $\beta=3.80$, chir. | $0.728716(56)$
$32\\!\times\\!32^{3}$, $\beta=3.80$, chir. | $0.728546(81)$
$24\\!\times\\!24^{3}$, $\beta=3.685$, symm. | $0.71209(13)$
$32\\!\times\\!32^{3}$, $\beta=3.80$, symm. | $0.738044(53)$
$56\\!\times\\!56^{3}$, $\beta=4.10$, symm. | $0.777900(17)$

The values for $Z_{\rm V}$ in tab. 3 were extracted from plateau ranges like
the one shown in fig. 5. Due to the simple expression for $Z_{\rm V}$ in eq.
(25), which only consists of a ratio of two correlation functions and yields a
good plateau quality in most cases, it can be determined to high precision. In
the present, preliminary results, plateaus have been taken from $x_{0}=T/4$ up
to $x_{0}=3T/4$.

##  7 Conclusions and outlook

We have seen that the Ward identity method in combination with Schrödinger
functional boundary conditions can be readily implemented to determine $c_{\rm
A}$ for the stabilised Wilson-Clover fermions. As we work in larger volumes,
some adaptions in the choice of the plateau range and wavefunctions turn out
to be advantageous to reach the desired precision goal. To finally arrive at a
proper interpolation over the whole range of interesting lattice spacings,
measurements on further ensembles are still to be added.

When it comes to renormalisation, a precise extraction of $Z_{\rm V}$ for the
vector current is straightforward, while the achievable accuracy for $Z_{\rm
A}$ remains to be seen. A promising route to enhance the final precision of
$Z_{\rm A}$ is to pursue the computational strategy within the $\chi$SF
advocated and successfully applied in [6, 14].

## References

  * [1] M. Lüscher, S. Sint, R. Sommer, and P. Weisz, Chiral symmetry and O(a) improvement in lattice QCD, Nucl. Phys. B 478 (1996) 365, [hep-lat/9605038].
  * [2] A. Francis, P. Fritzsch, M. Lüscher, and A. Rago, Master-field simulations of O(a)-improved lattice QCD: Algorithms, stability and exactness, Comput. Phys. Commun. 255 (2019) 107355, [1911.04533].
  * [3] J. Bulava and S. Schaefer, Improvement of $N_{f}$ = 3 lattice QCD with Wilson fermions and tree-level improved gauge action, Nucl. Phys. B 874 (2013) 188, [1304.7093].
  * [4] J. Bulava, M. D. Morte, J. Heitger, and C. Wittemeier, Non-perturbative improvement of the axial current in $N_{f}=3$ lattice QCD with Wilson fermions and tree-level improved gauge action, Nucl. Phys. B 896 (2015) 555, [1502.04999].
  * [5] J. Bulava, M. D. Morte, J. Heitger, and C. Wittemeier, Non-perturbative renormalization of the axial current in $N_{f}=3$ lattice QCD with Wilson fermions and tree-level improved gauge action, Phys. Rev. D 93 (2016) 114513, [1604.05827].
  * [6] M. Dalla Brida, T. Korzec, S. Sint, and P. Vilaseca, High precision renormalization of the flavour non-singlet Noether currents in lattice QCD with Wilson quarks, Eur. Phys. J. C 79 (2019) 23, [1808.09236].
  * [7] A. S. Francis, F. Cuteri, P. Fritzsch, G. Pederiva, A. Rago, A. Schindler, A. Walker-Loud, and S. Zafeiropoulos, Properties, ensembles and hadron spectra with Stabilised Wilson Fermions, PoS LATTICE2021 (2022) 118, [2201.03874].
  * [8] M. Della Morte, R. Hoffmann, and R. Sommer, Non-perturbative improvement of the axial current for dynamical Wilson fermions, JHEP 03 (2005) 029, [hep-lat/0503003].
  * [9] M. Della Morte, R. Hoffmann, F. Knechtli, R. Sommer, and U. Wolff, Non-perturbative renormalization of the axial current with dynamical Wilson fermions, JHEP 07 (2005) 007, [hep-lat/0505026].
  * [10] M. Lüscher, S. Sint, R. Sommer, P. Weisz, and U. Wolff, Non-perturbative O(a) improvement of lattice QCD, Nucl. Phys. B 491 (1996) 323, [hep-lat/9609035].
  * [11] U. Wolff, Monte Carlo errors with less errors, Comput. Phys. Commun. 156 (2003) 143, [hep-lat/0306017].
  * [12] F. Joswig, S. Kuberski, J. T. Kuhlmann, and J. Neuendorf, pyerrors: a python framework for error analysis of Monte Carlo data, 2209.14371.
  * [13] J. Heitger and F. Joswig, The renormalised $\mathrm{O}(a)$ improved vector current in three-flavour lattice QCD with Wilson quarks, Eur. Phys. J. C 81 (2021) 254, [2010.09539].
  * [14] S. Sint, The chirally rotated Schrödinger functional with Wilson fermions and automatic O(a) improvement, Nucl. Phys. B 847 (2010) 491, [1008.4857].