The temporal component of the potential and the divergent pieces of the radial potential have been given above in (<ref>) and (<ref>). One can clearly see by comparison with (<ref>) the difference between the BS and NU gauges. It should be noted that these two gauges do not only differ by $C$, but also by terms arising from the fact that the two solution spaces are actually different, as can be seen e.g. by comparing (<ref>) and (<ref>).

In the absence of boundary sources, the constant piece of the radial potential in NU gauge is given in (<ref>). Since this expression is quite involved, let us study it in the flat case and with $\partial_uq_{AB}=0=\delta\sqrt{q}$. For comparison with the BS gauge, we should furthermore rewrite the potential explicitly in terms of $C^{AB}_\tf$. This can be done using in particular (<ref>). Denoting $N_{AB}=\partial_uC^\tf_{AB}$, we then find with the simplifying assumptions that

\begin{equation}\label{NU r0 potential simple}
\setlength\fboxsep{10pt}
\boxed{
\begin{aligned}
\Theta^r_0\big|_\text{NU}
&=\f{1}{4}\sqrt{q}\left(2N_{AB}\delta C_\tf^{AB}-\delta q^{AB}\Big(C^\tf_{AB}R[q]-2D_AD^CC^\tf_{BC}+D_A\partial_BC\Big)+\delta CR[q]\right)\cr
\end{aligned}}
\end{equation}

We have omitted the total derivatives and variations since they do not play a role. Already in this simplified setup, where we have removed the cosmological constant, the boundary sources, the Weyl factor, and the time-dependency of the leading metric, we see that the trace $C$ appears in two places in the potential. It modifies the canonical pair conjugated to $q^{AB}$, and also appears as conjugated to the two-dimensional Ricci scalar $R[q]$. We therefore conclude that the finite piece of the potential in NU gauge differs from that in BS gauge by trace terms. These terms however drop if we further restrict the setup to that of the original BMS analysis, with a fixed two-sphere metric $q^\circ_{AB}$ such that $\delta q^\circ_{AB}=0=\delta R[q^\circ]$. Any relaxation of these historical boundary conditions will give rise to a difference between the BS and NU gauges.

The fact that the BS and NU potentials differ is a surprise, and it suggests that the two gauges may actually be physically inequivalent. This does not contradict the equivalence statements made e.g. in <cit.>, because these references set $C=0$. We already see on the simplified setup leading to (<ref>) that with $C\neq0$ there is an extra canonical pair in the phase space. The general formula (<ref>) shows that for a more general setup (e.g. with $\Lambda\neq0$ and a time-dependent Weyl factor) the BS and NU potentials will differ also by terms not involving $C$. Once again, this is because the solution space is different as we have seen above. As we are about to see in the next section, it turns out that the presence of $C$ in NU gauge is intimately related to the existence of an extra radial asymptotic symmetry. Ultimately, one should of course compute whether $C$ and/or the extra symmetry parameter appear in the charges for these asymptotic symmetries. We defer this heavy task to future work, and present for the time being the structure of the symmetries.

§.§ Vector fields and symmetries

In order to find the asymptotic Killing vectors in NU gauge, we simply need to study how the NU gauge condition (<ref>) constrains the expansion (<ref>) of the radial component found in the partial Bondi gauge. Preserving the NU gauge condition requires $\partial_r(\pounds_\xi g_{ur})=0$, where in NU gauge (<ref>) becomes
This in turn fixes the radial part of the vector field to satisfy

where $h(u,x^A)$ is an integration “constant”, which we have already shifted here in an $r$-independent manner in order to obtain a vector field forming an algebra instead of an algebroid. Upon further integrating (<ref>), we see that the radial component $\xi^r$ will contain two integration functions, $h(u,x^A)$ and a new $k(u,x^A)$, so that $\xi$ has in total 5 free functions $(f,Y^A,h,k)$ of $(u,x^A)$. To obtain the asymptotic expansion of the vector field we first use (<ref>) in NU gauge, and then integrate (<ref>) to get

ξ^r=r h+k+12e^2β_0(4∂_Aβ_0+D_A)∂^Af
where a log term has dropped because of (<ref>), and where we have once again used the freedom in redefining the integration constant in order to shift $k$. This expansion of the radial vector field should be compared with its BS analogue (<ref>). With our choices of integration constants, we see that these vectors conveniently differ at leading order only by the presence of the additional symmetry parameter $k$ in NU gauge. At lower orders the components $\xi^r$ in BS and NU gauge actually also differ, in a way which we do not make explicit here but which can easily be computed.

The transformation laws can be found by plugging the form of the radial vector field in NU gauge into the general formulas of section <ref>. Because at leading order the vector field is the same in BS and NU gauge up to the new symmetry parameter $k$, most of the transformation laws are also the same up to the action of $k$. More precisely, here in NU gauge we get

which agrees at leading order with the BS counterpart (<ref>). Note that here this equation is computing the transformation of the expansion (<ref>), while in BS gauge the transformation is exact because of the determinant condition. We see that here $h$ plays also the role of a Weyl rescaling, and that it is allowed to appear in (<ref>) precisely because we have considered the differential version (<ref>) of the NU gauge. The presence of the Weyl symmetry is therefore related to the inclusion of the boundary source $\beta_0$ in NU gauge, while it is related to the inclusion of a fluctuating boundary volume $\sqrt{q}$ in BS gauge. In either gauge, the Weyl symmetry is allowed because we are using a differential gauge condition rather than an algebraic one.  The algebra formed by the vector field in NU gauge, which contains the free functions $(f,Y^A,h,k)$, is the time-dependent algebra $\big(\text{diff}(\I^+)\loplus\mathbb{R}_h\big)\loplus\mathbb{R}_k$.

The transformation laws for the fields $\beta_0$, $U^A_0$, and $q_{AB}$ are the same as in BS gauge and given again by (<ref>), (<ref>), and (<ref>). The influence of the new symmetry generator $k$ is to act on subleading components of the angular metric, and we find

while for the trace we find

We therefore see that (<ref>) has acquired a shift in $k$ with respect to the BS analogue (<ref>).
The transformation laws (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>), are the generalization to four-dimensional gravity of the three-dimensional results (2.12) in <cit.>. Note that (<ref>) should be read with care since in BS gauge $C_{AB}$ is trace-free while this is not the case in NU gauge. What we mean by this equation is simply that the transformation of $C_{AB}$ in NU gauge is given by the same formula (<ref>) as the transformation in BS gauge, although now with a trace-full $C_{AB}$. In any case, this is consistent with plugging $\xi^r_{+1}$ and $\xi^r_0$ in (<ref>).

Equation (<ref>) shows that the new symmetry parameter $k$ appearing in the NU asymptotic Killing vector is related to the presence of the non-vanishing trace $C$ in the said gauge. This trace $C$ can be set to zero by a choice of origin for the affine parameter (which in NU gauge is the radial coordinate $r$), but the extra symmetry parameter at order $r^0$ in $\xi^r$ can shift this value. Said the other way around, one can argue that the presence of this symmetry allows to fix $C=0$. This is used explicitly e.g. in <cit.>. But the question remains whether this symmetry is gauge or physical, i.e. associated to a (non-)vanishing charge. In the three-dimensional analogue of the generalized NU gauge studied here, there is indeed a non-vanishing charge associated with these radial translations, and they furthermore form an Heisenberg algebra with the Weyl charges <cit.>. In order to study the (possible) charge associated with $k$ in the present four-dimensional context, one should first investigate the simplified solution space leading to the potential (<ref>) (where $C$ survives). The study of the Weyl sector is made much more complicated by the fact that in NU gauge it requires the presence of the boundary source $\beta_0$, which as we have shown in appendix <ref> renders the analysis of the symplectic potential extremely heavy. One could however simply use the Iyer–Wald or Barnich–Brandt formula (which can be applied without knowledge of the symplectic potential) to get a first idea of what the structure of the charges will be, but we expect, based on results from the three-dimensional theory, that extracting the proper charges will require the use of a corner ambiguity and also a possible change of slicing <cit.>. All these aspects will be investigated in forthcoming work.

§ PERSPECTIVES

In this article we have studied the structure of asymptotically locally (A)dS and flat gravity in Bondi coordinates, using the partial Bondi gauge (<ref>). This partial gauge, built with only three gauge fixing conditions and with no determinant condition, enables to treat together the Bondi–Sachs (BS) and Newman–Unti (NU) gauges, which can each be reached at any later stage of a calculation in the partial gauge by further imposing a specific fourth condition. So far, the Newman–Unti solution space has been studied in the literature using the Newman–Penrose formalism <cit.>. Here we have presented an analysis of the solution space which is by design similar to the analysis in BS gauge. This enables to clearly contrast the two gauge choices.

We have recalled why the NU gauge features an extra mode, which is the trace of $C_{AB}$, while this trace is set to vanish in BS gauge by the determinant condition. The NU gauge also admits an extra asymptotic symmetry, which is parametrized by a function at order $\O(r^0)$ in the radial part of the asymptotic Killing vector. Together with the supertranslations, the superrotations, and the Weyl transformations, which are also present in BS gauge, this extra component in NU gauge leads to the asymptotic symmetry algebra

The extra radial translation $k$ present in NU gauge acts precisely by shifting the trace $C$. The fate of this trace mode and that of the radial symmetry are therefore related. Indeed, one should now compute the asymptotic charges in order to see whether the radial symmetry (used to fix the value of $C$ in <cit.>) is pure gauge or physical (as in the three-dimensional case <cit.>). We have chosen to postpone this investigation to future work, but the key ingredients for this calculation have been laid out here.

In this study of the partial Bondi gauge we have also included logarithmic terms in the solution space by working with the polyhomogeneous expansion (<ref>). We have seen that, consistently, all these logarithmic terms are forced to vanish when $\Lambda\neq0$ when completing the partial gauge to BS or NU gauge and going on-shell. When $\Lambda=0$ however, the logarithmic terms are involved in evolution equations and also modify the Weyl scalars. Moreover, we have kept in the construction of the solution space a free boundary metric with time-dependent fields $\beta_0$, $U^A_0$, and $q_{AB}$. The solution space constructed here is therefore an extension of that in <cit.> in two respects: the relaxation to the partial Bondi gauge, and the inclusion of the logarithmic terms. The prospects are however similar, namely to also study sourced holography and understand radiation in (A)dS. We have seen how, in this endeavor, the Newman–Penrose formalism using the Weyl scalars and spin coefficients proves to be very useful. This has enabled us in particular to recover the (A)dS evolution equation (<ref>) for the mass, and to identify the quantity $\S^{AB}$ sourcing this equation and conjugated to $q_{AB}$ in the symplectic structure (<ref>), as the radiation tensor shifted by the covariant spin-2. This use of the NP formalism bypasses the symmetry arguments used in <cit.> to identify the covariant functionals and derive the evolution equations.

There are many interesting prospects to be explored in future work, and which rely on some of the extensions of the asymptotic structure and the solution space presented here. A non-exhaustive list is as follows:

* Radial translation charges.
An immediate and important follow-up is to compute the asymptotic charges and the asymptotic charge algebra in the partial Bondi gauge. The question is that of the fate of the new generator $k$ in NU gauge, or equivalently the term of order $\O(r^0)$ in the partial Bondi gauge expansion (<ref>). In three dimensions, we have shown in <cit.> that this radial translation is a physical symmetry with a non-vanishing charge, and a choice of corner term and integrable slicing was necessary in order to obtain integrable charges. In the four-dimensional case, one can expect that similar subtleties will arise, and that there will be a preferred genuine slicing in which the charges become integrable in the absence of radiation. Although, as we have seen, there is no flux-balance law associated with the trace mode $C$, it would be interesting to investigate whether and how this changes in the presence of matter <cit.> and in modified theories of gravity <cit.>.

* Subleading Weyl scalars and higher spin charges.
In the flat case there is an infinite tower of functionals whose time evolution is constrained. These are the subleading components of the angular metric, which appear in “covariant” form in the subleading terms $\Psi^{n>5}_0$ in $\Psi_0$. We have seen here how $E^\tf_{AB}$ appears in covariant form $\E_{AB}$ in $\Psi_0^5$. These subleading terms are related to the higher spin charges <cit.> and to the Newman–Penrose conserved charges <cit.>. We have explained here how these subleading terms can have mixed helicities. It would be interesting for future work to analyse these subleading terms in details (not necessarily in the partial Bondi gauge), in order to understand their relationship with the higher spin symmetries and to relate them with the operators defined in the context of celestial holography.

* Logarithmic terms.
Another aspect of the present work which should be explored further concerns the polyhomogeneous expansion and the role of the logarithmic terms. While in the asymptotically locally (A)dS case the logarithmic terms are necessarily vanishing on-shell, in the flat case they are subject to evolution equations which we have derived. The question is now how these logarithmic terms affect the symplectic renormalisation of the potential and the construction of the charges. We have also seen how the logarithmic terms relax the standard fall-offs of Weyl scalars, and one should now study how they modify the subleading terms in the Weyl scalars and the Newman–Penrose charges <cit.>. Moreover, it would be interesting to understand whether the logarithmic terms in the solution space are related to the log soft theorem <cit.>. Finally, let us mention an example where logarithmic modes play an important role. This is the case of three-dimensional chiral gravity <cit.>, where it has been shown that there is a new mode, the “logarithmic primary”, which violates the usual Brown–Henneaux boundary conditions and requires a larger class of boundary conditions <cit.>. While the setup is quite different from the present analysis, this shows nonetheless that logarithmic modes can arise in certain contexts. One should also extend the present construction to include matter in order to see whether the logarithmic terms can lead to new flux-balance laws and memory effects.

* Radiation in (A)dS.
The notion of radiation for spacetimes with $\Lambda\neq0$ is poorly understood <cit.>. Although here we have not directly discussed notions pertaining to gravitational waves and quadrupole formulas in (A)dS, we have seen how certain well-defined quantities in asymptotically flat spacetimes generalize to the case $\Lambda\neq0$ and $\partial_uq_{AB}\neq0$ (the question of whether these remain the relevant quantities is however open). For example, we have proposed a notion of generalized news in (<ref>), and also identified the role of the unconstrained spin-2 in the quantity $\S^{AB}$ sourcing the mass loss and entering the symplectic structure. We hope to clarify in future work the role of these quantities when studying the memory effects in (A)dS.

* Robinson–Trautman solutions.
Our extended solution space, which includes in particular arbitrary time-dependent sources for the boundary metric, is adapted to the study of the algebraically special Robinson–Trautman spacetimes <cit.>. These spacetimes are interesting in their own right in order to study sourced holography and memory effects, but could also provide a simplified setup in which to study the role of the trace mode $C$ and the radial symmetry in NU gauge.

§ ACKNOWLEDGMENTS

It is our pleasure to thank Miguel Campiglia, Adrien Fiorucci, Laurent Freidel, Daniel Grumiller, Yannick Herfray, Alok Laddha, Pujian Mao, Blagoje Oblak, Marios Petropoulos, Daniele Pranzetti, Ana-Maria Raclariu, Romain Ruzziconi, Ali Seraj, and Simone Speziale for numerous discussions. We also warmly thank Miguel Campiglia and the Universidad de la República for their hospitality during the Montevideo Workshop on Celestial Symmetries, where part of this work was completed. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities.

§ APPENDICES

§.§ Notations and conventions

We gather here some notations and conventions used throughout the main text, as well as some useful identities.

* All the angular indices $A,B,\dots$ are lowered and raised with the leading order angular metric $q_{AB}$. We therefore have $C^{AB}=q^{AC}q^{BD}C_{BD}$, and $U_A=q_{AB}U^B$.

* The one and only exception to this rule is with the full angular metric $\gamma_{AB}$, whose true inverse such that $\gamma^{AC}\gamma_{BC}=\delta^A_B$ is denoted $\gamma^{AB}$.

* It is sometimes more compact to lower and raise indices with other two-dimensional tensors, in particular in the study of the symplectic potential in appendix <ref>, so in these cases we write

* The contraction of indices is denoted with square brackets, and in particular we have

as well as objects such as $[CC]=C^{AB}C_{AB}$.

* The traces in the metric $q_{AB}$ are denoted by index-free objects, e.g. $C=q^{AB}C_{AB}$.

* Trace-free components are denoted by $C^\tf_{AB}=C_{AB}-\f{1}{2}Cq_{AB}$, and it is easy to see that we therefore have $[C^\tf D^\tf]=[CD^\tf]=[C^\tf D]$.

* The covariant derivative associated with $\gamma_{AB}$ is $\D_A$ and satisfies $\gamma_{AB}\D_CU^B=\D_C\tU_A$.

* The covariant derivative associated with $q_{AB}$ is $D_A$ and satisfies $q_{AB}D_CU^B=D_CU_A$.

* For the expansions in $r$ we use the big-$\O$ notation, so that solving an equation $\text{E}=\O(r^{-n})$ means that the terms in $r^{-n+1}$ and above have been killed. By abuse of notation we also use $\O(r^{-n})$ even for the polyhomogeneous expansions containing logarithmic terms.

* We define symmetrization as $P_{(A}Q_{B)}=P_AQ_B+P_BQ_A$, i.e. without factor $1/2$.

A few useful identities used throughout the main text are
C_ABδq^AB+C^ABδq_AB   =0,

q_ABδC^AB+q^ABδC_AB   =2δC,

C_(ACD^C_B)   =([CD]-CD)q_AB+DC_AB+CD_AB,

C_ACC^C_B   =12([CC]-C^2)q_AB+CC_AB,

C_ACC^CDC_DB   =12C([CC]-C^2)q_AB+12([CC]+C^2)C_AB,

12δ[CC]   =(C^2-[CC])δln√(q)+C^AB(δC_AB-Cδq_AB),

(δln√(q)+δ)∂_uln√(q)   =1√(q)δ∂_u√(q),

=∂_uC^_ABδC_^AB-12C^_ABδ(C∂_uq^AB)   +(12q_AB∂_u+∂_uq_AB)(CδC_^AB)   +12(12C^2∂_uq_ABδq^AB+∂_uCδC+CδC∂_uln√(q)-C∂_uCδln√(q)).
We give other useful identities involving covariant derivatives in appendix <ref> below.

§.§ Christoffel coefficients

We gather here the components of the Christoffel connection for the metric (<ref>) in partial Bondi gauge. Noting that here the repeated angular indices in $\Gamma^A_{\mu A}$ are summed, and using the notation $\tU_A=\gamma_{AB}U^B$, we have

Γ^u_μr   =0,

Γ^A_rr   =0,

Γ^u_uu   =12e^-2β∂_r(e^2βVr+γ_ABU^AU^B)+2∂_uβ,

Γ^u_Au   =-12e^-2β∂_r_A+∂_Aβ,

Γ^u_AB   =12e^-2β∂_rγ_AB,

Γ^r_uu   =12e^-2β((Vr∂_r+U^A∂_A)g_uu+∂_uγ_ABU^AU^B)+12r(2V∂_uβ-∂_uV),

Γ^r_ur   =-12e^-2β(∂_r(e^2βVr)+12γ_AB∂_r(U^AU^B))-U^A∂_Aβ,

Γ^r_rr   =2∂_rβ,

Γ^r_uA   =-12e^-2β(Vr∂_r_A+U^B∂_uγ_AB+U^B∂_B_A+_B∂_AU^B))-12r∂_AV,

Γ^r_Ar   =12e^-2βγ_AB∂_rU^B+∂_Aβ,

Γ^r_AB   =12e^-2β((∂_u+Vr∂_r)γ_AB+_A_B+_B_A),

Γ^A_uu   =12U^A(e^-2β∂_rg_uu+4∂_uβ)-12γ^AB(∂_Bg_uu+2∂_u_B),

Γ^A_uB   =-12e^-2βU^A∂_r_B+12γ^AC(∂_C_B-∂_B_C+∂_uγ_BC)+U^A∂_Bβ,

Γ^A_uA   =-14e^-2β(∂_rγ_ABU^AU^B+∂_r[U])+∂_uln√(γ)+U^A∂_Aβ,

Γ^A_rB   =12γ^AC∂_rγ_BC,

Γ^A_rA   =∂_rln√(γ),

Γ^A_ur   =12γ^AB(∂_Be^2β-∂_r_B),

Γ^C_AB   =Γ^C_AB[γ]+12e^-2βU^C∂_rγ_AB.

§.§ Spin coefficients and frame identities

Using the vectors $\ell$ and $n$ defined in (<ref>) and the null frames (<ref>), the Ricci rotation coefficients are defined by
Below we give their explicit expression in the partial Bondi gauge, as well as their on-shell expansion using the expansion (<ref>) for the frames. We find
α   =12(γ_124-γ_344)=12(V2rΓ^u_rA-Γ^r_rA-m̂̅̂^B∇_Am̂_B)m̂̅̂^A=12rD_Am̅_0^A+Ø(r^-2),

β   =12(γ_123-γ_343)=12(V2rΓ^u_rA-Γ^r_rA-m̂̅̂^B∇_Am̂_B)m̂^A=-12rD_Am_0^A+Ø(r^-2),

γ   =12(γ_122-γ_342)=14r^2e^-2β(r∂_rV-V)-12n^μm̂̅̂^A∇_μm̂_A=Λ6r+γ_0+γ_1r+Ø(r^-2),

ϵ   =12(γ_121-γ_341)=18∂_rγ_AB(m̂̅̂^Am̂̅̂^B-m̂^Am̂^B)-∂_rβ   =12(γ_121-γ_341) =18r^2C_AB(m̅_0^Am̅_0^B-m_0^Am_0^B)+Ø(r^-3),

π   =γ_421=-(12e^-2βγ_AB∂_rU^B+∂_Aβ)m̂̅̂^A=1r^2(e^-2β_0U_A^2+C_AB∂^Bβ_0)m̅_0^A+Ø(r^-3),

ν   =γ_422=(e^-2βV4r(2∂_Aβ-e^-2βγ_AB∂_rU^B)-Γ^r_μAn^μ)m̂̅̂^A   =γ_422 =Λ3r∂_Aβ_0m̅_0^A+(12e^-2β_0∂_AV_+2m̅_0^A+Λ3∂_Aβ_0m̅_1^A)+Ø(r^-1),

μ   =γ_423=-12e^-2β(∂_uln√(γ)+V2r∂_rln√(γ)+_AU^A)=-Λ6r-Λ24C+μ_1r+Ø(r^-2),

λ   =γ_424=-12e^-2β((∂_u+V2r∂_r)γ_AB+2_A_B)m̂̅̂^Am̂̅̂^B

=γ_424 =-12(Λ6C^_AB-e^-2β_0𝔅^Λ_AB)m̅_0^Am̅_0^B+Ø(r^-1),

κ   =γ_131=0,

τ   =γ_132=(∂_Aβ-12e^-2βγ_AB∂_rU^B)m̂^A=2r∂_Aβ_0m^A_0+Ø(r^-2),

σ   =γ_133=12∂_rγ_ABm̂^Am̂^B=-12r^2C_ABm_0^Am_0^B+Ø(r^-3),

ρ   =γ_134=12∂_rln√(γ)=1r-C4r^2+Ø(r^-3),
=14e^-2β_0(V_+2+m_0^A∂_A(U^0_Bm̅_0^B)-m̅_0^A∂_A(U^0_Bm_0^B)-m̅^0_BU^A_0D_Am_0^B)   +18e^-2β_0(m̅_0^Am̅_0^B-m_0^Am_0^B)(𝔅^Λ_AB+2D_AU_B^0-Λ6e^2β_0C^_AB),

=Λ96(4D-[CC])+Λ24(14CC^_AB-𝔏_AB)(m^A_0m^B_0-m̅^A_0m̅^B_0)   +18e^-2β_0(∂_uq_AB-𝔅^Λ_AB)(m^A_0m^B_1-m̅^A_0m̅^B_1)+(D_Am̅^A_0m^B_0-D_Am^A_0m̅^B_0)∂_Bβ_0   +18e^-2β_0(N_AB-12C𝔅^Λ_AB)(m^A_0m^B_0-m̅^A_0m̅^B_0),

The frames (<ref>) used in the NP formalism satisfy important identities. First we have the time derivative

∂_um_0^A   =(14∂_uq_AB(m̅_0^Am̅_0^B-m_0^Am_0^B)-12∂_uln√(q))m_0^A-12(∂_uq_BCm_0^Bm_0^C)m̅_0^A,
which can in turn be used in (<ref>) to compute the time evolution $\partial_um^A_1$. Then, we have other identities involving the time evolution, which are

m̅_0^Am̅^B_0(∂_uq_AC)q^CD(∂_uq_DB)   =23m̅_0^Am̅^B_0(ΛC_AB-3𝔅^Λ_AB)∂_uln√(q),

m̅_0^Am̅^B_0(∂^2_uq_AB)   =m̅_0^Am̅^B_0(Λ3(∂_u+∂_uln√(q))C_AB-(∂_u+∂_uln√(q))𝔅^Λ_AB).
Finally we have the identities
C_ABm^A_0m^B_0   =C_AB^m^A_0m^B_0=iq_ACϵ^CDC^_DBm^A_0m^B_0,

m_0^BD_Bm_0^A   =m_0^AD_Bm_0^B,

m_0^BD_Bm̅_0^A   =-m̅_0^AD_Bm_0^B,

(m_0^A∂_A+2β_1)(V_Am^A_0)   =(D_AV_B)m^A_0m^B_0,

(m^A_0∂_A-4β_1)(T_ABm̅^A_0m̅^B_0)   =(D^BT^_AB)m̅^A_0,

2σ_2V_Am̅^A_0   =-C^_ABV^Am^B_0,

(V_Am̅^A_0)(T_ABm^A_0m^B_0)   =T^_ABV^Am^B_0,

2 Re((m_0^A∂_A-2β_1)(V_Am̅^A_0))   =D_AV^A,

4 Re(σ_2T_ABm̅^A_0m̅^B_0)   =-C^_ABT^AB=-C_ABT_^AB,
which hold for any vector $V_A$ and any symmetric tensor $T_{AB}$, and which can easily be complex-conjugated to obtain mirror relations.

To compute the above expansions of the spin coefficients we have also used the expansion of the normal $n^\mu=(n^u,n^r,n^A)$, which is

n^u   =e^-2β_0+Ø(r^-2),

=r^2Λ6+r2e^-2β_0V_+2+12e^-2β_0V_+1+Λ96(4D-[CC])   +12re^-2β_0(V_0+196([CC]-4D)(6∂_uln√(q)+6D_AU^A_0-e^2β_0ΛC))+Ø(r^-2),

n^A   =e^-2β_0U^A_0+2r∂^Aβ_0+Ø(r^-2).

§.§ Expansion of the inverse angular metric

Here we explain how to obtain analytically the asymptotic expansion of the determinant and the inverse of the angular metric. We initially assume that there is no logarithmic term in the expansion of the angular metric, but the construction can straightforwardly be applied to take such terms into account.

Let us denote all the matrices appearing in (<ref>) by $q\equiv q_{AB}$, $C\equiv C_{AB}$, and so on. We also denote by $q^{-1}\equiv q^{AB}$ the inverse to $q$, by $C^*\equiv C^{AB}=q^{-1}Cq^{-1}$ the tensor with indices raised by $q_{AB}$, and by a dot the matrix product $C^*\!.D^*=C^*qD^*$ in $q$. We then use the fact that, for arbitrary matrices $A$ and $B$ such that $A$ and $A+B$ are non-singular, we have
Applying this formula to $\gamma$ (and discarding once again for simplicity the logarithmic terms) with $A=r^2q$ and $B=rC+D+\O(r^{-1})$ enables to extract the coefficients in
In particular, this gives
γ^*_4   =-D^*+C^*.C^*,

γ^*_5   =-E^*+C^*.D^*+D^*.C^*-C^*.C^*.C^*,

γ^*_6   =-F^*+D^*.D^*+C^*.E^*+E^*.C^*-C^*.C^*.D^*-C^*.D^*.C^*-D^*.C^*.C^*+C^*.C^*.C^*.C^* .
This can then be rewritten using the fact that for $2\times2$ matrices we have
AB+BA   =(A)B+(B)A-(A+B)𝕀+(A)𝕀+(B)𝕀,

(A+B)   =A+B+(A)(A^-1B),
where $\tr(\cdot)$ is the usual matrix trace. Introducing now the trace with respect to $q$ defined as $\tr_q(C)\equiv\tr(q^{-1}C)=q^{AB}C_{AB}$ and using the notation $[CC]\equiv\tr(C^*C)=C^{AB}C_{AB}$ we find
(q^-1C)   =12(_q(C)^2-[CC]),

(q^-1C)(C^-1D)   =(q^-1D)(D^-1C)=[CD]-_q(C)_q(D),
which we can use to obtain
C^*.C^*   =12([CC]-_q(C)^2)q^-1+_q(C)C^*,

C^*.D^*+D^*.C^*   =([CD]-_q(C)_q(D))q^-1+_q(D)C^*+_q(C)D^*,

C^*.C^*.C^*   =12(_q(C)([CC]-_q(C)^2)q^-1+([CC]+_q(C)^2)C^*).
Combining these formulas and switching back to the notation of the main text finally yields

γ_4^AB   =-D^AB+CC^AB_+12[CC]q^AB,

The subleading components $\gamma_4^{AB}$ and $\gamma_5^{AB}$ are unfortunately necessary because they appear in the symplectic potential and when solving respectively for $\beta_2$ and $\beta_3$ in section <ref>.

To obtain the determinant, we proceed with the same notations as above and use formula (<ref>) to expand
=r^4(q)+r^2(C)+(D)+(Er+Ø(r^-2))   +(q) [q^-1(r^3C+r^2D+rE+D+Fr+Ø(r^-2))]   +(C)[C^-1(rD+E+Dr+Ø(r^-2))]   +(D)[D^-1(Er+Ø(r^-2))].
To go back to the notations of the main text we then note that
and similarly with the other tensors appearing in the angular metric.

Combining all these results, we can for example extract the next order in the partial Bondi gauge solution (<ref>), which is used in appendix <ref>.

At the end of the day, we can also reintroduce the logarithmic terms in order to expand the determinant and the inverse of the polyhomogeneous metric (<ref>). In this case we find that the angular metric and the volume element have an expansion of the form

=r^2+C2r+14(2D-[CC^])+12r(E+14C[CC^]-[CD^]+E_llnr)   +14r^2(132(3C^4-4[CC]^2-4C^2[CC])-34C^2D+12D[CC]+12D^2+C[CD]-[DD].   .+CE-2[CE]+2F+(CE_l-2[CE_l]+2F_l)lnr+2F_l^2(lnr)^2)+Ø(r^-3).
All of these expressions can of course be obtained with a computer, although there can be a lot of guess work involved as more and more combinations of the tensors appear as one goes to lower order in the radial expansion.

§.§ Expansion of the angular connection

In this appendix we assume once again for simplicity and clarity that there are no logarithmic terms in the expansion of the angular metric. The Levi–Civita connection associated with the angular metric $\gamma_{AB}$ then admits an expansion of the form
Γ^C_AB[γ]   =Γ^C_AB[q]+Γ^C_1ABr+Γ^C_2ABr^2+Γ^C_4ABr^3+Ø(r^-4),
Γ^C_1AB   =12(D_AC^C_B+D_BC^C_A-D^CC_AB),

=12q^CD(∂_AD_DB+∂_BD_AD-∂_DD_AB)-12C^CD(∂_AC_DB+∂_BC_AD-∂_DC_AB)   +γ^4_MNq^MCΓ^N_AB[q],

=12q^CD(∂_AE_DB+∂_BE_AD-∂_DE_AB)-12C^CD(∂_AD_DB+∂_BD_AD-∂_DD_AB)   +12γ_5^CD(∂_AC_DB+∂_BC_AD-∂_DC_AB)+γ^5_MNq^MCΓ^N_AB[q],
with $\gamma^4_{AB}$ and $\gamma^5_{AB}$ given above in (<ref>). When the determinant condition is satisfied we have at all orders the properties

or in terms of the above expansion $\Gamma^B_{nAB}\stackrel{\text{\tiny{BS}}}{=}0$, $\forall n$.

We now give a few identities involving combinations of covariant derivatives. In the computation of (<ref>) we use identities of the form
D_(AP̅_B)-C_(ACD_B)P^C-2Γ^C_1ABP_C   =P^CD_CC_AB,

D_AP̅^A-C^ABD_AP_B-Γ^A_1P_A   =Γ^B_1ABP^A=12P^A∂_AC,
which hold for any vector $P^A$ with $\bar{P}^A=C^{AB}P_B$, and where we have denoted $q^{AB}\Gamma^C_{1AB}=\Gamma^C_1$. With the expansion (<ref>) we have
For any symmetric trace-free tensor $T_{AB}$ we have

Finally we note that

which also holds when replacing e.g. $\delta\to\partial_u$ or $\delta\to\partial_A$.

§.§ Subleading solution space

Here we report part of the structure of the solution space at subleading order. This is necessary in order to compute the evolution equations in the case $\Lambda\neq0$ and also to understand how the logarithmic branches disappear in the case $\Lambda\neq0$. Since we compute the evolution equations in section <ref> with the boundary sources turned off, i.e. $\beta_0=0=U^A_0$, we also do so in this appendix.

The solution for $\beta$ is actually independent of the boundary sources $\beta_0$ and $U^A_0$. Going further in the expansion of (<ref>), we find terms in $r^{-4}\big(\beta_{l^24}(\ln r)^2+\beta_{l4}\ln r+\beta_4\big)$, with
β_l^24   =-38F_l^2,

β_l4   =332CE_l+932[CE^_l]-38F_l+14F_l^2,

32 β_4
=-12F+4F_l+9[EC^]+3CE-CE_l-154[E_lC^]+4[𝔏D]+D^2   +18D(3C^2-22[CC])-5C[𝔏C]-332C^2[CC]+58[CC]^2.
The expansion (<ref>) then continues with terms in $r^{-4}\big(U^A_{l^24}(\ln r)^2+U^A_{l4}\ln r+U^A_4\big)$, with
U^A_l^24   =0,

U^A_l4   =-34(12C U^A_l3+C^A_BU^B_l3-D_B(E_l^)^AB-13∂^A E_l),

U^A_4   =(lengthy).
Finally, we find that the expansion (<ref>) continues with $r^{-1}\big(V_{l^21}(\ln r)^2+V_{l1}\ln r+V_1\big)$, where
V_l^21   =5Λ12F_l^2,

V_1   =(lengthy).
The expressions for $U^A_4$ and $V_1$ have been used for the computation of the NP (A)dS mass loss equation, but are too lengthy to be displayed here.

We now give the real parts of the subleading terms in the Weyl scalar $\Psi_4$ in the absence of boundary sources. Off-shell of the $(AB)$ Einstein equations we have

=-12D^A_A-34C+Λ24C_AB^AB   +Λ48(5[CE^_l]+[C𝔏])-12D_AD_B𝔏^AB+14[N𝔏]-18[C𝔏]∂_uln√(q).

On-shell of the constraints and in the case (A)dS we therefore get $\text{Re}\big(\Psi_2^{l4}\big)=0$ and (<ref>), which is what we have used in the main text.

§.§ Alternative resolution of $\boldsymbol{\mathrm{E}_{ur}=0}$

The field equation $\mathrm{E}_{ur}=0$ can rewritten using its trace to obtain the equivalent form $R_{ur}=\Lambda g_{ur}$. At the difference with $\gamma^{AB}R_{AB}=2\Lambda$, which is first order, this equation is second order in $V$ since we have
=(∂_r-∂_rln√(γ))Γ^r_ur+∂_AΓ^A_ur+∂_u∂_r(ln√(γ)+2β)   +14(∂_uγ_AB)(∂_rγ^AB)+12γ^ABΓ^C_AC[γ](∂_Be^2β-∂_r_B),
which contains $\partial_r\Gamma^r_{ur}$. The expansion which solves this equation is of the form
which differs from (<ref>) by the presence of an additional $r\ln r$ term. The rewritten equation sets
R_ur-Λg_ur=Ø(r^-1)   ⇒ V_+3=Λ3e^2β_0,

R_ur-Λg_ur=Ø(r^-2)   ⇒ V_+2=Λ6e^2β_0C-∂_uln√(q)-D_AU^A_0,

R_ur-Λg_ur=Ø(r^-3)   ⇒ V_+l1=-12[C𝔅^Λ],

R_ur-Λg_ur=Ø(r^-4)   ⇒ V_l0=[D𝔅^Λ]-14C[C𝔅^Λ]+Λ3e^2β_0E_l,
where the logarithmic terms involve $\mathfrak{B}^\Lambda_{AB}$ defined in (<ref>). As expected, we see that with this method of resolution the logarithmic terms involve the constraints whose vanishing is implied by the $(AB)$ Einstein equations. The construction of the solution space is therefore consistent if we follow the route presented in this appendix or that of section <ref>.

§.§ Computation of $\boldsymbol{\Theta^r}$

Here we present the details of the calculation of the radial part $\Theta^r$ of the symplectic potential in the partial Bondi gauge, assuming there are no logarithmic terms in the angular metric and the solution space. We first rewrite the radial component of the potential in the form
= √(γ)(δ(Γ^u_uu+Γ^A_uA-Γ^r_ur)+U^Aδ(Γ^u_Au+Γ^B_AB-Γ^r_Ar)+VrδΓ^A_rA+e^2βγ^ABδΓ^r_AB)   √(γ)(δ𝔸+U^Aδ𝔹_A+ℂ+𝔻).
We then study each term separately using the off-shell expansion of the connection components (<ref>). The only simplification which we use is to set $\beta_1=0$ since at the end of the day this on-shell condition is common to both the BS and NU gauges. We also assume that the logarithmic branches have been removed so that e.g. $U^A_{l3}$ does not appear in the expansion of $U^A$.

The first term in the above rewriting is
𝔸_+1   =2V_+3,

𝔸_0   =-e^-2β_0[U_0U_1]+2[U_0∂β_0]+V_+2+∂_u(ln√(q)+2β_0),

=-e^-2β_0(3([U_0U_3]+[U_1U_2])+[_1U_1]+2[_0U_2]+[_0U_1]-2β_2[U_0U_1])   +2([U_0∂β_2]+[U_2∂β_0])-2M-4β_2V_+2-6β_3V_+3+14∂_u(2D-[CC]+8β_2).The second term is
𝔹_0A   =e^-2β_0U^1_A+∂_Aln√(q),

𝔹_1A   =e^-2β_0(2U^2_A+^1_A)+12∂_AC,

𝔹_2A   =e^-2β_0(3U^3_A+2^2_A+^1_A-2β_2U^1_A)+14∂_A(2D-[CC]).
The third term is
2ℂ_0   =-V_+3δC,

2ℂ_1   =-V_+2δC+V_+3δ([CC]-2D),

2ℂ_2   =-V_+1δC+V_+2δ([CC]-2D)+12V_+3δ(C^3-6E-3C[CC]+6[CD]).
Finally, the fourth and most intricate term is

𝔻_+1   =q^ABδ(q_ABV_+3)-4δβ_0V_+3,

𝔻_0   =q^ABδ{q_ABV_+2-12C_ABV_+3+12∂_uq_AB+D_AU^0_B}   +δβ_0(CV_+3-4V_+2-∂_ulnq-2D_AU^A_0)+V_+3δC,

=q^ABδ{q_ABV_+1-12C_ABV_+2+12∂_uC_AB+D_AU^1_B+C_ACD_BU^C_0+12U^C_0D_CC_AB}   -C^ABδ{12∂_uq_AB+D_AU^0_B}   -δβ_0(4V_+1-CV_+2+([CC]-2D)V_+3+∂_uC+2D_AU^A_1+U^A_0∂_AC)   +V_+2δC+12δV_+3([CC]-2D)+V_+3(γ_4^ABδq_AB-12C^ABδC_AB-4δβ_2),

=q^ABδ{2Mq_AB+12C_ABV_+1+12∂_uD_AB-12E_ABV_+3}   +q^ABδ{D_A(U^2_B+^1_B+^0_B)-Γ^C_1AB(U^1_C+^0_C)-Γ^C_2ABU^0_C}   -C^ABδ{q_ABV_+1+12C_ABV_+2+12∂_uC_AB+D_A(U^1_B+^0_B)-Γ^C_1ABU^0_C}   +γ_4^ABδ{q_ABV_+2+12C_ABV_+3+12∂_uq_AB+D_AU^0_B}+γ_5^ABδ(q_ABV_+3)   -δβ_0(8M-CV_+1-2DV_+2-EV_+3+q^AB∂_uD_AB+[CC]V_+2-C^AB∂_uC_AB)   -δβ_0(γ_4^AB(C_ABV_+3+∂_uq_AB+2D_AU^0_B)+2γ_5^ABq_ABV_+3)   -2δβ_0(D_A(U^A_2+^A_0-C^AB^0_B)+12(U^A_1+^A_0)∂_AC+(C^ABΓ^C_1AB-Γ^C_2)U^0_C)   +δβ_2(CV_+3-4V_+2-∂_ulnq)-4δβ_3V_+3.
We then collect all these pieces and use the expansion (<ref>) of the boundary volume to find

=δ𝔸_2+𝔻_2+ℂ_2+[U_0δ𝔹_2]+[U_1δ𝔹_1]+[U_2δ𝔹_0]+12Cθ^r_1+18(4D-C^2-2[CC])θ^r_2   +18(4E+3C[CC]-4[CD]-2CD-12C^3)θ^r_3.
Here we clearly see in each $\theta^r_n$ the terms involving $\theta^r_{m>n}$ which come from the relaxation of the determinant condition, or in other words the fact that the co-dimension 2 volume has a non-trivial expansion. It is immediate to see that these terms drop in the BS gauge.

These expressions are so far valid off-shell. The only assumptions we have made about the expansion is $\beta_1=0$ and the absence of logarithmic terms. It is now convenient to use the on-shell expressions derived in the partial Bondi gauge in section <ref> in order to separate these components of the potential into sources pieces (S) which vanish when the boundary sources $(\beta_0,U^A_0)$ are set to zero, and fixed pieces (F) which remain even in the absence of boundary sources. For compactness and convenience we will however adopt a mixed notation where we still keep e.g. $U^A_1$ without its on-shell value since it is a pure source term. We also keep $\beta_2$ and $\beta_3$ unspecified in order to later be able to reach the BS or NU gauge, and sometimes keep writing $V^+_n$ instead of  their explicit on-shell values.

We now decompose each coefficient in the expansion of the potential as $\theta^r_n=\o{F}\theta^r_n+\o{S}\theta^r_n$, keeping in mind that we are now using the solution space of section <ref>. For the terms which remain in the absence of boundary sources we find
øFθ^r_3   =V_+3δlnq,

øFθ^r_2   =(2CV_+3-∂_ulnq-∂_u)δln√(q)-12δq^AB∂_uq_AB+12V_+3(q_ABδC^AB+2δC),

=-12q^ABδ{q_AB[e^2β_0(R[q]+Λ2(34[CC]-C^23-D))+12(∂_uln√(q)+2∂_u)C]}   +14q^ABδ{C_AB(∂_ulnq-CV_+3)+2∂_uC_AB}-12C^ABδ∂_uq_AB   +14(CV_+3-∂_ulnq+2∂_u)δC+V_+3(γ_4^ABδq_AB-12C^ABδC_AB+14δ[CC])   +12CøFθ^r_2+18(4D-C^2-2[CC])øFθ^r_3,

=2(δlnq+δ)M+12q^ABδ{∂_uD_AB-E_ABV_+3+D_A(e^2β_0(∂_BC-D^CC_BC))}   +14q^ABδC_AB{e^2β_0(R[q]+Λ2(34[CC]-C^23-D))+12(∂_uln√(q)+2∂_u)C}   -14C^ABδ{C_AB(CV_+3-∂_ulnq)+2∂_uC_AB}+12(δCV_+1-CδV_+1)   +12γ_4^AB(δ{q_AB(CV_+3-∂_ulnq)}+δC_ABV_+3+δ∂_uq_AB)+γ_5^ABδq_ABV_+3   +14∂_uδ(2D-[CC]+8β_2)+2δ{β_2(∂_ulnq-CV_+3)}+δβ_2(∂_ulnq-CV_+3)   +14(CV_+3-∂_ulnq)δ([CC]-2D)+14V_+3δ(C^3-6E-3C[CC]+6[CD])   +12e^2β_0(∂^AC-D_BC^AB)∂_Aδln√(q)+12CøFθ^r_1+18(4D-C^2-2[CC])øFθ^r_2   +18(4E+3C[CC]-4[CD]-2CD-12C^3)øFθ^r_3-10δβ_3V_+3.

For the terms which drop when turning off the sources we find
øSθ^r_3   =2δV_+3,

=2CδV_+3-δ{e^-2β_0[U_0U_1]-2[U_0∂β_0]-2∂_uβ_0+3D_AU^A_0}+[U_0δ𝔹_0]   +q^ABδ(D_AU^0_B)-(D_AU^A_0)δlnq-δβ_0(CV_+3-∂_ulnq-2D_AU^A_0),

=-12q^ABδ{q_AB[e^2β_0(4D_A∂^Aβ_0+8(∂^Aβ_0)(∂_Aβ_0))+12(D_AU^A_0+2U^A_0∂_A)C]}   +q^ABδ{D_AU^1_B+C_ACD_BU^C_0+12U^C_0D_CC_AB}-C^ABδ(D_AU^0_B)   +18δV_+3(3[CC]-4D)-12(D_AU^A_0)δC+12q^ABδ(C_ABD_CU^C_0)   -δβ_0(4V_+1-CV_+2+([CC]-2D)V_+3+∂_uC+2D_AU^A_1+U^A_0∂_AC)   -δ{e^-2β_0([U_1U_1]+2[U_0U_2]+[_0U_1])}+2δ[U_1∂β_0]   +[U_0δ𝔹_1]+[U_1δ𝔹_0]+12CøSθ^r_2+18(4D-C^2-2[CC])øSθ^r_3,

=-q^ABδ{D_A(e^2β_0C_BC∂^Cβ_0)-D_A(^1_B+^0_B)+Γ^C_1AB(U^1_C+^0_C)+Γ^C_2ABU^0_C}   +14q^ABδC_AB{e^2β_0(4D_A∂^Aβ_0+8(∂^Aβ_0)(∂_Aβ_0))+12(D_AU^A_0+2U^A_0∂_A)C}   +C^ABδ{12C_ABD_CU^C_0-D_A(U^1_B+^0_B)+Γ^C_1ABU^0_C}   +γ_4^AB(12C_ABδV_+3-δ(q_ABD_CU^C_0)+δ(D_AU^0_B))+γ_5^ABq_ABδV_+3   -δβ_0(8M-CV_+1-2DV_+2-EV_+3+q^AB∂_uD_AB+[CC]V_+2-C^AB∂_uC_AB)   -δβ_0(γ_4^AB(C_ABV_+3+∂_uq_AB+2D_AU^0_B)+2γ_5^ABq_ABV_+3)   -2δβ_0(D_A(U^A_2+^A_0-C^AB^0_B)+12(U^A_1+^A_0)∂_AC+(C^ABΓ^C_1AB-Γ^C_2)U^0_C)   -δ{e^-2β_0(3([U_0U_3]+[U_1U_2])+[_1U_1]+2[_0U_2]+[_0U_1]-2β_2[U_0U_1])}   +2δ([U_0∂β_2]+[U_2∂β_0])-4δ(β_2D_AU^A_0)+2δβ_2D_AU^A_0-12D_AU^A_0δ([CC]-2D)-6β_3δV_+3   +[U_0δ𝔹_2]+[U_1δ𝔹_1]+U^A_2δ(e^-2β_0U^1_A)-e^2β_0C^AB∂_Aβ_0∂_Bδln√(q)+12CøSθ^r_1   +18(4D-C^2-2[CC])øSθ^r_2+18(4E+3C[CC]-4[CD]-2CD-12C^3)øSθ^r_3.

From these expressions it is then possible to reduce the potential to the BS or NU gauges. This reduction is done by imposing either of the following  two conditions:
(BS)   ←  C=0,  D=12[CC],E=[C],

(NU)   ←   C=0,  D=14[CC],E=23([C]+18C(C^2-2[CC])),
where in NU gauge the conditions on $D$ and $E$ are equivalent to $\beta_2=0$ and $\beta_3=0$ respectively.

In the main text we have focused on the potential $\o{F}\theta_0$ in the absence of boundary sources, i.e. with $\beta_0=0=U^A_0$. We have furthermore split the potential between a part which survives in BS gauge, and additional terms which are only present in NU gauge (involving essentially $C$). The divergent pieces are reported in (<ref>) and (<ref>). The constant piece $\o{F}\theta_0$ is however still very lengthy so we report it here. In BS gauge it is given by

=2(δlnq+δ)M-12q^ABδ(Λ_AB+D_AD^CC_BC)   +14q^ABδC_AB(R[q]+5Λ8[CC])+12δC^AB∂_uC_AB-18∂_u(C^ABδC_AB)   +12(32C^AB∂_uC_AB-D_AC^AB∂_B)δln√(q)   +18([CC]δln√(q)-3C^ABδC_AB)∂_uln√(q).
In NU gauge we find instead

=øFθ^r_0|_BS   +(Λ12Cγ_4^AB-Λ12(C^2+[CC])C^AB+12CC^AB∂_uln√(q))δq_AB   +12(37Λ72C^3-11Λ18C[CC]+C(C∂_u-R[q]-C^AB∂_uq_AB)+∂^AC∂_A-14[CC]∂_uln√(q))δln√(q)   -(Λ18C+132(∂_ulnq+2∂_u))δ[CC]   +(Λ6C^2-2Λ9[CC]-14∂_uC+78C∂_uln√(q))δC   +(Λ6CC^AB-Λ6C^2q^AB+7Λ96[CC]q^AB+18q^AB(∂_uln√(q)+2∂_u)C)δC_AB   +12q^ABδ(D_A∂_BC)+18q^ABδ∂_u(CC^_AB)-12δ(CC^AB∂_uq_AB)   -14δ(CR[q])+18Cδ((2∂_u+∂_uln√(q))C)-18∂_u(C^2δln√(q)-CC^ABδq_AB)   +18√(q)([CC]+C^2)δ∂_u√(q)-Λ9√(q)δ(√(q) [C])-116q^ABδ(∂_uq_AB[CC]+q_AB∂_u[CC])   +14√(q)C(∂_u(√(q) q^ABδC_AB)-∂_uδ(√(q) C)-∂_u(δ√(q) C)).
This expression can surely be simplified further, but we refrain from doing so here. One should just note that the additional terms which differentiate the BS and NU potential are not only proportional to $C$. The reason for this is two-fold. First, we see by comparing (<ref>) and (<ref>) that the on-shell expression of the subleading components of the angular metric differ in BS and NU gauge. Second, we should recall that in BS gauge $\beta_{n>0}\neq0$ while in NU gauge $\beta_{n>0}=0$. Since the NU potential is written in (<ref>) as the sum of the BS potential and extra terms, these terms need to remove the non-vanishing $\beta_{n>0}\neq0$ which are present in the BS potential (in order to ensure that for the NU potential we have $\beta_{n>0}=0$).

§.§ Off-shell $\boldsymbol{\mathrm{E}^\tf_{AB}\big|_{\O(r^{-2})}}$

We report here the off-shell expression for $\mathrm{E}^\tf_{AB}\big|_{\O(r^{-2})}
^{\beta_0=0=U_0}$ . This is the expression from which the equations (<ref>), (<ref>), and (<ref>) have been derived in the main text (upon further imposing $\Lambda=0$ and going on-shell of $\mathfrak{B}^0_{AB}=0$). The length of this expression also illustrates the clear advantage of the NP formalism over the use of the tensorial Einstein equations. We have

2Λ3F_l^2AB^(lnr)^2-EOM(E_lAB)lnr   +Λ{13F_l^2AB^-56F_lAB^+23F^_AB-14C_(AC(E_l^)^C_B)+16CE^_AB   +Λ{+132(43D[CC^]+56C^2[CC]-16C^4-[CC]^2)q_AB   +Λ{+118C^_(AC𝔏^CDC^_DB)-19[CC^]𝔏_AB+18C[C𝔏]q_AB-16[D𝔏]q_AB}
-23{∂_u_AB-12(D_(A_B))^-32C^_AB|_Λ=0+32q_ACϵ^CDC^_DB   -23{+12(_AB-3E^_l AB-38C𝔏_AB-94C^_ABD+158𝔏_(AC(C^)^C_B))∂_uln√(q)   -23{+12(∂_uC+32R[q])𝔏_AB+18C∂_u𝔏_AB+34(C^)_(A^C∂_u𝔏_B)C   -23{+916D(C^)^C_(A∂_uq_CB)+716𝔏_(ACC^CD∂_uq_DB)-118(C^)_(A^C(∂_uq_CD)𝔏^D_B)   -23{+316C_(AC^𝔏^CD𝔅^0_DB)-532C𝔏_(A^C𝔅_CB)^0+14[𝔏𝔅^0]C^_AB   -23{+2132C_(ACC_^CDC_DE^(𝔅^0)^E_B)-332[CC]C^_(AC(𝔅^0)^C_B)   -23{-4164C C_(ACC^CD𝔅^0_DB)+C[CC^]𝔅^0_AB+716C^2C^_(AC(𝔅^0)^C_B)-316C[C𝔅^0]C^_AB}
-13((2∂_u+∂_uln√(q))E_lAB^+12(D_(AU_B)^l3)^)   +(12E^C_l(A-4E^C_(A)𝔅^0_B)C+5(E-14E_l)𝔅^0_AB,
with $\M\big|_{\Lambda=0}=M+\f{1}{16}\big(\partial_u+\partial_u\ln\sqrt{q}\big)\big(4D-[CC]\big)$ and

+E^_l(AC𝔅^0C_B)+34E_l𝔅^Λ_AB+12[E_l𝔅^Λ]q_AB   +Λ3(5F^_l^2AB-2F^_lAB-12CE_lAB^+18([CC]-C^2)𝔏_AB+18(C_(AC𝔏^CDC_DB))^).We note that the first term on the second line is always vanishing on-shell. This is because we have $\mathfrak{B}^0_{AB}=0$ in the flat case and $E^\tf_{lAB}=0$ in (A)dS.


# Automated Real-Space Lattice Extraction for Atomic Force Microscopy Images

Marco Corrias  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
<EMAIL_ADDRESS>University of Vienna, Vienna Doctoral School in
Physics, Boltzmanngasse 5, 1090 Vienna, Austria  Lorenzo Papa  University of
Vienna, Faculty of Physics and Center for Computational Materials Science,
Kolingasse 14-16, A-1090 Vienna, Austria  Igor Sokolović  Institute of Applied
Physics, TU Wien, Wiedner Hauptstrasse 8-10/134, A 1040 Vienna, Austria
Viktor Birschitzky  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
University of Vienna, Vienna Doctoral School in Physics, Boltzmanngasse 5,
1090 Vienna, Austria  Alexander Gorfer  University of Vienna, Faculty of
Physics and Center for Computational Materials Science, Kolingasse 14-16,
A-1090 Vienna, Austria  Martin Setvin  Department of Surface and Plasma
Science, Faculty of Mathematics and Physics, Charles University, V
Holešovičkách 2, 180 00 Prague 8, Czech Republic  Institute of Applied
Physics, TU Wien, Wiedner Hauptstrasse 8-10/134, A 1040 Vienna, Austria
Michael Schmid  Institute of Applied Physics, TU Wien, Wiedner Hauptstrasse
8-10/134, A 1040 Vienna, Austria  Ulrike Diebold  Institute of Applied
Physics, TU Wien, Wiedner Hauptstrasse 8-10/134, A 1040 Vienna, Austria
Michele Reticcioli  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
Cesare Franchini  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
Dipartimento di Fisica e Astronomia, Università di Bologna, 40127 Bologna,
Italy

###### Abstract

Abstract. Analyzing atomically resolved images is a time-consuming process
requiring solid experience and substantial human intervention. In addition,
the acquired images contain a large amount of information such as crystal
structure, presence and distribution of defects, and formation of domains,
which need to be resolved to understand a material’s surface structure.
Therefore, machine learning techniques have been applied in scanning probe and
electron microscopies during the last years, aiming for automatized and
efficient image analysis. This work introduces a free and open source tool
(AiSurf: Automated Identification of Surface Images) developed to inspect
atomically resolved images via Scale-Invariant Feature Transform (SIFT) and
Clustering Algorithms (CA). AiSurf extracts primitive lattice vectors, unit
cells, and structural distortions from the original image, with no pre-
assumption on the lattice and minimal user intervention. The method is applied
to various atomically resolved non-contact atomic force microscopy (AFM)
images of selected surfaces with different levels of complexity: anatase
TiO2(101), oxygen deficient rutile TiO2(110) with and without CO adsorbates,
SrTiO3(001) with Sr vacancies and graphene with C vacancies. The code delivers
excellent results and has proved to be robust against atom misclassification
and noise, thereby facilitating the interpretation scanning probe microscopy
images.

##  I Introduction

In recent years giant leaps have been made in scanning probe microscopy,
particularly in atomic force microscopy (AFM) Giessibl (2003); Morita  _et
al._ (2015); Pavliček and Gross (2017). Atomically sharp tips, often
functionalized with simple molecules (e.g., CO-terminated tips) Scheuerer  _et
al._ (2019); Kempkes  _et al._ (2019); Fournier  _et al._ (2011); Wagner  _et
al._ (2012) allow not only for the structural identification of material
surfaces at the atomic level but also for precise manipulation of single atoms
and molecules Gross  _et al._ (2009); Alldritt  _et al._ (2020).
Improved imaging techniques, aided by artificial intelligence (AI) Krull  _et
al._ (2020); Kalinin  _et al._ (2021a, 2022), lead to a vast abundance of
atomically resolved images. Acquisition rates have indeed reached a point
where the analysis of individual images by humans, albeit computer-assisted,
is proving to be a new bottleneck in the advancement of related surface
science fields.

Machine learning (ML) can represent a viable alternative to accelerate the
processing of these data. Supervised learning approaches based on deep neural
networks (DNN) have been applied to a variety of different tasks Choudhary
_et al._ (2022), ranging from tip functionalization Rashidi and Wolkow (2018);
Gordon  _et al._ (2020) and lattice recognition Ziatdinov  _et al._ (2018);
Kalinin  _et al._ (2021b); Ziletti  _et al._ (2018); Leitherer  _et al._
(2021), all the way up to automated molecular structure discovery Alldritt
_et al._ (2020); Carracedo-Cosme  _et al._ (2021). Despite these promising
results, the applicability of DNN for efficient interpretation of experimental
images is limited by the available datasets Kalinin  _et al._ (2016). While
electron microscopy (EM) Aversa  _et al._ (2018); Schwenker  _et al._ (2020);
Ede (2020) and scanning tunneling microscopy (STM) Choudhary  _et al._ (2019)
image datasets have been recently published, we are not aware of publicly
available dataset of atomically resolved experimental AFM images for surfaces
(a database on simulated AFM images for molecular identification has been
recently collected Carracedo-Cosme  _et al._ (2022)). Augmenting an image
dataset using measured and simulated data are sub-optimal options for
practical purposes: (i) although AI optimization methods have been proposed
Kalinin  _et al._ (2021a), it is complicated to measure high-quality images
and define a reliable reproducibility protocol; (ii) simulating atomically
resolved images using computational surface science methods requires the
precise knowledge of the structural model at the atomic scale, which might be
accessible for bulk-terminated surfaces but highly complex in the case of
surface structural reconstructions; the vast diversity of defects that can be
present in a surface complicates this task even further; (iii) last but not
least, supervised machine learning also heavily relies on the network
architecture, which forbids non-expert users to play with different
architectures to achieve better accuracies. Open-source packages like
pycroscopy Somnath  _et al._ (2019) and AtomAI Ziatdinov and Kalinin (2021)
can overcome this initial difficulty, although they remain dependent on a
dataset. These reasons hinder the adoption of supervised machine learning in
the field of AFM microscopy.

In the absence of large datasets, unsupervised machine learning methods
constitute a convenient alternative to data-hungry supervised ML techniques
and slow manual analysis of individual images. Several successful attempts
based on unsupervised algorithms have been made. Atomap Nord  _et al._ (2017)
analyzes scanning transmission electron microscopy (STEM) images by using 2D
Gaussian fitting, differentiating atomic columns by their different
brightness; different approaches based on clustering and principal component
analysis (PCA) were also proposed Belianinov  _et al._ (2015); Somnath  _et
al._ (2018). In a recent work focused on the automated identification of local
structures in atomically resolved images, Laanait et al. Laanait  _et al._
(2016) have shown that the well-established computer vision algorithms Scale-
Invariant Feature Transform Lowe (2004) (SIFT) combined with Clustering
Algorithms (CA) Gan  _et al._ (2007) are capable of recognizing and labeling
atoms according to their local environment. An advantage of this approach is
that it does not involve any pre-assumptions on the underlying periodicity of
the lattice. Reciprocal space analysis can be used to analyze periodic
structures in images and, eventually, remove image artifacts. However, some
prior knowledge of the system and the intervention by an expert user are
required. Restricting the analysis to real space is expected to provide better
stability under such features, thereby expanding the degree of analysis.

In this work, we propose an open-source tool, AiSurf (”Automated
identification of Surface images”), developed to inspect and classify
crystalline 2D phases in atomically-resolved images via Scale Invariant
Feature Transform and Clustering Algorithms. With no pre-assumption on the
lattice symmetry and only minimal human intervention, AiSurf can be applied
for the analysis of regular structures, as well as for images showing complex
structural correlations involving the formation of structural domains,
structural and chemical defects, and recognition of adsorbates. To assess the
performance and transferability of the code, we analyze a variety of surface
structures characterized by different features, specifically: (i) Defect-free
anatase TiO2(101); (ii) SrO-terminated SrTiO3(001) with Sr vacancies; (iii)
Rutile TiO2(110) with O vacancies; (iv) CO adsorbates on rutile TiO2(110); (v)
Simulated graphene with C vacancies. Additionally, to further test the
capabilities of the algorithm against experimental noise, we show the results
obtained for an SrTiO3(001) image affected by noise and other artifacts. All
AFM images have been acquired via AFM experiments, apart from the graphene
image, which has been simulated using the probe particle model Hapala  _et
al._ (2014, 2016).

The algorithm and computational protocol are presented in the next section.
The experimental details on the acquisition of the AFM images and their
automated analysis are presented and discussed in Sec. III. Data, code and
documentation are available as indicated in Sec. V.

##  II Method

Figure 1: Summary of the workflow of AiSurf using an atomic force microscopy
image of the anatase TiO2(101) surface. a) AFM image with SIFT keypoints
colored according to their cluster label. b) Distance vectors within the
(blue) reference cluster. Colored circles show the result of a clustering done
in this figure; each cluster represents a lattice vector candidate: $\vec{a}$
and $\vec{b}$ have been chosen. c) Distance vectors connecting keypoints. The
colors of the vectors represent the deviation from the predicted lattice
vectors; the top-right inset shows the distribution of keypoints in the unit
cell, with keypoints colored according to their cluster; the bottom-left inset
shows the predicted unit cell with its sublattice positions. We point out that
what is interpreted as Ti atom is just a minimum, since Ti atoms are too far
down to be reached by the tip.

This section describes the workflow of AiSurf, graphically schematized in Fig.
1.

First, we apply SIFT (as implemented in the library ’opencv’ Bradski (2000))
to an input image to extract atomic features (”keypoints”) as bright and dark
blobs, in analogy with the work of Laanait et al. Laanait  _et al._ (2016). We
specify that ”atomic features” may also refer to adsorbates on the surface or,
in general, to minima or maxima present in the image. SIFT detects keypoints
using the difference-of-Gaussians (DoG) function with sub-pixel resolution.
Although, for optimal detection, the input parameters have to be adjusted
(depending on the scale and contrast of the features in the image), AiSurf
default values can handle the images satisfactorily. Every parameter is
described in the documentation provided in Sec. V. Ideally, all visible
features should be identified while no keypoints should be assigned to
background noise. To further enhance the quality of keypoints, we exclude a
thin area on the border of the image, and we ignore keypoints showing a
significant deviation in size from the median value. This constraint should
ensure that only atomically-sized features remain.

The keypoints detected by SIFT are characterized by a feature description
based on image gradients that are brightness, scale, and rotational invariant.
We break rotational invariance to distinguish keypoints belonging to different
sublattices if they show the same local environment but different
orientations. This can be the case of atoms of the same species but in
different sublattices (e.g., graphene lattice, see Fig. 2), or atoms next to a
point defect (i.e. vacancy) but in different orientations.

The keypoints are then clustered based on their SIFT descriptors by using
agglomerative clustering taken from the scikit-learn library Pedregosa  _et
al._ (2011). This method is deterministic and only uses the number of clusters
as an input parameter. Density-based clusterings are challenging to use here
because they are typically parameterized by a characteristic distance, which
is difficult to estimate for abstract descriptors like ours, which are a
collection of local gradients. The optimal number of clusters is chosen by
calculating the silhouette score Rousseeuw (1987) for each number of clusters
in a chosen interval and taking the one that maximizes it. A label is assigned
to each cluster.

As an example, we show the results obtained for a defect-free AFM image of the
anatase TiO2(101) surface. The experimental image overlaid with the clustered
keypoints is shown in Fig. 1a. The clustering labels allow us to select a
reference cluster that will be used to extract the lattice vectors. It should
contain only keypoints on the same sublattice since they have similar local
environments and, thus, similar SIFT descriptors. We point out that AiSurf is
currently not designed to distinguish different domains; for example,
translational or rotational domains Hommrich  _et al._ (2002); Lee  _et al._
(2006) will be recognized as a single one. The reference cluster is by default
the one with most keypoints, but users can change it at their preference. In
order to obtain the lattice vectors, the nearest neighbours are computed for
every keypoint of the reference cluster. These nearest neighbors are grouped
into sub-clusters, used to calculate the distance vectors from the reference
point. The centers of mass of these sub-clusters represent possible lattice
vector candidates. The two shortest, linear independent candidates are chosen
as lattice vectors $\vec{a}$ and $\vec{b}$. Finally, a Bravais lattice is
generated by linear combinations of these lattice vectors. Fig. 1b shows a
visual representation of this process. We note that the identification of
lattice vectors and Bravais lattice is based uniquely on the reference
cluster, and all other clusters are not verified to hold the same vectors:
this choice allows to exclude spurious effects arising from defects.

After determining the lattice vectors, it is possible to calculate the
sublattice positions by processing all detected keypoints again. For each
keypoint we compute its position within the newly-obtained unit cell; the
resulting distribution of positions within the unit cell is then clustered.
The center of mass of each cluster is calculated using the Bai-Breen method
Bai and Breen (2008), accounting for the periodic boundary conditions. The
resulting distribution of keypoints in the unit cell is visible in the top-
right inset of Fig. 1c. Up to this point, the keypoints were labeled according
to the clustering performed on the descriptors, which allowed for the
selection of the reference cluster for the lattice vectors identification.
Now, every keypoint is labeled according to its sublattice position.

The algorithm can also plot an average view of the unit cell by averaging the
area around each keypoint belonging to the same atomic species. The predicted
unit cell and its sublattice positions are then drawn. The bottom left of Fig.
1c shows the image extracted with this method.

##  III Results and discussion

###  III.1 Experimental and simulated AFM images

The images of anatase TiO2(101) in Fig. 1, the SrO termination of cleaved
SrTiO3(001) in Fig. 2b,c and rutile TiO2(110) in Fig. 3 have been obtained
using non-contact atomic force microscopy (nc-AFM) Morita  _et al._ (2015);
Giessibl (2019) in constant-height mode and ultra-high vacuum, whereas the nc-
AFM image of graphene was simulated. The signal displayed as grayscale image
is the frequency shift $\Delta$f. Details on the experimental and
computational setups are given below.

Fig. 2a shows a constant-height nc-AFM image simulated using the probe
particle model Hapala  _et al._ (2014, 2016) on free-standing graphene with C
vacancies in a regular pattern, and three N atoms substituting C in the sites
neighboring the vacancy. The simulation was performed using a CO-terminated
AFM tip Gross  _et al._ (2009) (oriented with the O atom towards the surface)
at a distance of 0.36 nm from the most-protruding atom in the surface slab,
using an oscillation amplitude of 200 pm. The CO tip, with $-0.05e^{-}$ charge
on the protruding O atom, was characterized by a spring constant of 0.5 N/m.
The simulation was performed over a $8\times 8$ nm2 surface area, with
$802\times 802$ pixels.

Figure 2b shows an unreconstructed SrO-terminated region on a cleaved
SrTiO3(001) surface Sokolović  _et al._ (2019, 2021). This termination
consists of Sr and O atoms, organized in a (1$\times$1) square pattern, with a
characteristic concentration of 0.14$\pm$0.02 monolayers of point defects in
the form of Sr vacancies. The image was acquired in close vicinity to the
surface using an O-terminated tip Sokolović  _et al._ (2020), such that Sr and
O atoms are detected in the attractive (dark) and repulsive regime (bright),
respectively; Sr vacancies are imaged as bright, cross-shaped features. The
gradually decaying contrast from the bottom to the top of the image is due to
the over-compensation of vertical drift caused by creep of the piezo scanner,
moving the AFM sensor away from the surface. The image was obtained with no
application of tip-sample bias voltage, with an oscillation amplitude of 100
pm, over a $5.4\times 5.4$ nm2 surface region with $400\times 400$ pixels.
Imaging was performed in ultra-high vacuum with base pressure below
1$\times$10-11 mbar, at a sample temperature of 5 K.

Fig. 2c shows the same surface as Fig. 2b. In addition to the intrinsic Sr
vacancies imaged in the form of cross-shaped point defects, this surface also
hosts protruding defects (in a concentration of less than 0.5%, displayed as
dark regions) that typically appear after a surface is exposed to the residual
gas in ultra-high vacuum for several days. The contrast differences in the
four image regions are due to different tip-sample distances throughout the
image since the nc-AFM tip was manually retracted/approached to the surface to
avoid contact with the protruding defects. The image was obtained with an
oscillation amplitude of 300 pm, without applying a bias voltage, over a
16.5$\times$16.5 nm2 region with 400$\times$400 pixels. Imaging was performed
under the same conditions of Fig. 2b.

The nc-AFM image of clean rutile TiO2(110) shown in Fig. 3a was acquired
during a previous study Reticcioli  _et al._ (2017). This surface consists of
rows of two-fold coordinated bridging O2c atoms and five-fold coordinated Ti5c
atoms, running along the $[001]$ and alternating along the $[1\bar{1}0]$
direction; the rows of O atoms are occasionally interrupted by single point
defects in the form of oxygen vacancies. Bright spots indicate O atoms
detected in the repulsive regime, whereas O vacancies are imaged as missing
spots (the Ti atoms are too distant from the tip and are not detected by nc-
AFM).

The constant-height nc-AFM image of a rutile TiO2(110) covered by CO molecules
shown in Fig. 3b was acquired during a previous study Reticcioli  _et al._
(2019). Imaging was performed with a CO-terminated tip that detects each
adsorbed CO molecule as a bright spot; bridging oxygen atoms are not resolved
because the adsorbed CO molecules protrude significantly more.

###  III.2 Automated analysis and discussion

Figure 2: Application of AiSurf to different atomic force microscopy images:
simulated graphene with C vacancies (a, d); SrTiO3(001) with Sr vacancies (b,
c, e, f). Panels a–c: Surfaces with detected keypoints; color coding indicates
different clusters. On graphene, C atoms appear bright (repulsive). On SrTiO3,
Sr vacancies are imaged as bright, cross-shaped features, Sr atoms are imaged
as dark spots, O atoms as bright spots. Panels d–f: Corresponding lattice
analysis for different sublattices (the color gradients indicates deviations
from the predicted perfect lattice); the bottom insets show the predicted unit
cell. The top-right panels in (d) and (e) are schematic views of the surface
structures.

The results obtained by AiSurf for graphene and SrTiO3(001) are collected in
Fig. 2. The parameters used for this analysis are included in the code
repository.

Three images were analyzed to test the capabilities of AiSurf under different
conditions. Fig. 2a shows a simulated image of graphene, used here to inspect
the capability of our algorithm in fully controlled conditions; Fig. 2b shows
a typical experimental AFM image of SrTiO3(001) with high contrast, which
represents a more realistic case; the experimental image in Fig. 2c shows
multiple artifacts, and strong contrast variations; it has been used to test
the capabilities of the algorithm under challenging conditions.

Figs. 2a–c show (part of) the keypoints detected by AiSurf, marked with
coloured circles (each cluster with a different color). These are the
keypoints left after the filtering process described in Sect. II. It can be
noticed how the SIFT algorithm can accurately detect the centers of the
features, even if they are only a few pixels in size, as the ones in Fig. 2a,
c. This is vital for the lattice recognition analysis since off-centered
keypoints lead to inaccurate predictions. This first clustering process may
not correctly detect all the different atomic species in the image but
guarantees that the largest clusters contain features of the same type, which
is needed for a successful analysis. Indeed, Figs. 2a–c show that atoms
surrounding defects are typically assigned to different clusters for the same
atomic species far from the defects; this apparent misclassification gets
corrected during the sublattice recognition, where keypoints are labelled
according to their sublattice position instead of the descriptor.

Figure 3: Automated analysis of experimental atomic force microscopy images of
rutile TiO2(110) with (b, d) and without (a, c) CO adsorbates. Panel a:
bridging oxygen atoms protruding from the rutile surface show rows of bright
spots interrupted by oxygen vacancies; the red circles indicates the automatic
detection of these oxygen atoms. The inset shows a sketch of the surface
structure: the protruding bridging oxygen atoms are labelled as ’O2c’, while
the surface Ti and O atoms are labelled as ’Ti5c’ and ’O3c’, respectively. AFM
image adapted from Ref. Reticcioli  _et al._ (2017). Panel c: deviation plot,
showing the deviations from the lattice vectors. The predicted unit cell is
shown in the inset. Panel b: adsorbed CO molecules are shown as bright spots;
CO on Ti5c atoms are marked with red circles; AFM image adapted from Ref.
Reticcioli  _et al._ (2019). Panel d: analysis of the CO-CO distance for every
CO on Ti5c atom detected in the experimental image.

In all cases, the identification of lattice vectors and unit cells works with
good accuracy. Since defects are neglected in the lattice vector extraction,
they pose no challenge to the algorithm. However, some atoms are not detected,
leading to gaps in the lattices shown in the bottom row images. In addition,
the breaking of rotational invariance leads to a separation of the carbon
sublattices (see Fig. 2a) and a distinction between keypoints next to a defect
(based on the direction the defect lies). The knowledge of the lattice
periodicity allows us to identify any deviation from a regular pattern. Fig.
2d–f show the vectors between keypoints belonging to the same sublattice, with
a color coding depicting their deviation from the predicted average lengths.
This distance-deviation plot helps detect displacements due to the presence of
defects and other sources. Fig. 2d shows slight displacements only around the
nitrogen impurities. The top-right inset depicts the ideal lattice and
highlights the two inequivalent C atoms (”C1” and ”C2”). The bottom-left inset
shows the average predicted unit cell; an orange and a blue circles represent
the two carbon atoms. The algorithm has recognized the green circle at the
center of the cell as a feature due to its strong contrast with the other
ones, but it is known that such feature is an atom for hexagonal lattices, not
honeycomb ones like this case. All the images, except for the insets with
white background, have been automatically generated by the algorithm.
Fig. 2e shows the Sr and O sublattices and their distance deviations. As
expected, deviations are present in the proximity of defects; ’Sublattice 2’
shows some misclassified points but evident distortions. The square unit cell
has been flawlessly predicted. Fig. 2f, used to test the algorithm’s
capabilities, shows overall positive results. Sublattice 1, defined by Sr
atoms, presents several unrecognized areas, especially in the low-contrast
regions. Distance deviations might not be reliable in this case. Sublattice 2
posed a lower challenge for its detection. The zoomed areas show that the
algorithm can detect features even when barely visible, proving its high
robustness against noise. The unit cell has been successfully predicted.

Fig. 3 shows the application of AiSurf to rutile TiO2(110) with and without CO
adsorbates Reticcioli  _et al._ (2019). Fig. 3a shows some detected bridging
oxygen atoms, highlighted in red. The corresponding deviation plot is shown on
Fig. 3c: distance deviations near oxygen vacancies are visible; away from
vacancy sites, no relevant deviations are present. This image is a good
example of how defect-induced deviations can be easily highlighted with
AiSurf. Such local structural distortions could play functional role in
chemical reactions between surface and adsorbates and might also facilitate
the formation and identification of metastable configurations Meier  _et al._
(2022). The inset in Fig. 3c shows the average unit cell. Other than the
central oxygen atom (labeled as ”O2c” in Fig. 3a inset), two other features
are present in the cell, interpreted as three-fold coordinated oxygen atoms
(”O3c”) given their position with respect to O2c.

Figs. 3b,d show the analysis of CO molecule adsorbed on rutile TiO2(110). CO
molecules adsorb both on five-fold coordinated Ti atoms and oxygen vacancies.
CO on Ti atoms are highlighted with red circles in Fig. 3b. The CO-CO
distances distribution on these sites is shown in panel (d). In the latter one
can easily visualize details barely visible from the original AFM image:
neighboring CO along $\vec{a}$ ([001] direction) tend to repel each other by
slightly tilting away from the surface normal, unlike second nearest neighbor
CO or CO belonging to different rows. This can be noticed by observing the
different shape of the two clusters at the center of panel (d), which are
differently distributed than the others, marking an absence of short-distance
neighboring CO. This result agrees with density functional theory calculations
Kunat  _et al._ (2009); Prates Ramalho  _et al._ (2017). Without computer
tools, such details can be overlooked in content-rich images like Fig. 3b.

While AiSurf is an aid for finding features in atomically resolved AFM images,
determination of unit cells, and deviations from a regular arrangement, the
interpretation is mostly left to the user. First, it has to be noted that the
SIFT algorithm is based on detection of local minima and maxima. A minimum or
maximum is not necessarily a physical feature, as exemplified by the hollow
sites in the hexagonal graphene rings if Fig. 2a,d or the minima in Fig. 1c
(blue dot in the bottom inset, no atom is detected in this position, Ti atom
is an interpretation). In the SrTiO3 case, the maxima between the dark Sr
sites (sublattice 1 in Fig. 2e,f) are at the positions of the O atoms, but
this is a mere coincidence. Areas with Sr vacancies show that the O sublattice
is not resolved; these areas appear with roughly constant brightness and do
not show the O atoms as maxima.

When analyzing displacements, it is important to be aware of the influence of
the tip on imaging. Especially tips with a rather wobbly termination such as
CO-functionalized tips can easily deform due to (electrostatic) forces on the
tip Hapala  _et al._ (2014); this can lead to apparent displacements of the
surface atoms Hapala  _et al._ (2016). Oxygen-terminated tips tend to be
stiffer, but even in the absence of tip deformation the minima or maxima in
$\Delta$f images may not exactly coincide with the atomic positions. In our
SrTiO3 example, the Sr vacancies are charged defects, thus they will distort
the electrostatic field in their vicinity. Since the negative frequency shift
above the Sr atoms is due to electrostatic attraction between the O-terminated
tip and the positive Sr atom, the field of a neighboring vacancy will
influence the position of the $\Delta$f minimum. Thus, the displacements
around the Sr vacancies (yellow and orange vectors in Fig. 2e) do not
necessarily reflect the atomic coordinates, and they may be also related to
the imaging process.

##  IV Summary and Conclusions

In this work, various atomically resolved non-contact AFM images with
different surface symmetries and degrees of complexity were analyzed in an
automated way using the open source tool AiSurf. At its present state, AiSurf
allows detecting the distribution of interatomic distances in distinct
sublattices and extracts primitive lattice vectors, unit cell, and distance
deviations from ideally symmetric lattices. We have shown results obtained on
simulated graphene with vacancies, as well as experimental images of anatase
and rutile TiO2 (the latter with and without CO adsorbates) and SrTiO3(001)
with Sr vacancies. The distance-deviation plot proved to be useful for the
detection of lattice distortions caused by oxygen vacancies, and the CO
tilting (see Fig. 3c, d respectively). This could be applied to the analysis
of any sources of bond deviation and similiar. The algorithm performed well
even on a SrTiO3(001) experimental image strongly affected by noise and
artifacts. Such robustness suggests a future implementation for real-time
analysis, parallel to the image acquisition process. For this reason, the
application on atomically resolved images other than AFM ones (e.g. scanning
tunneling microscopy) is believed to be possible and will be further tested.
In addition, future developments are already planned, such as the
identification and analysis of grain boundaries, domains and detection of
defects.

AiSurf is a user-friendly, documented unsupervised-ML tool that requires
minimal user intervention and no need to provide any image database.
Developing a robust tool for analyzing atomically resolved images is still an
ambitious goal for both supervised and unsupervised techniques. The most
adopted ones deliver successful results for limited types of images, which
usually present a good atomic contrast and/or are acquired with a precise
technique, such as the one for which they were designed. We aim to design a
protocol that relies on little or no prior knowledge of the physical system,
hoping to achieve a wide degree of flexibility. To conclude, we believe that
providing intuitive and open access tools will help material scientists
accelerate time-consuming tasks like image selection and analysis, and assist
the detection of elusive features in atomically-resolved images.

##  V Data availability statement

The AiSurf tool along with the necessary input data will be soon publicly
available.

## Acknowledgments

This work was supported by the Austrian Science Fund (FWF) project Super
(Grant No. P 32148-N36) and the SFB F81 project TACO. The computational
results have been achieved using the Vienna Scientific Cluster (VSC). MSe
acknowledges the support from the Czech Science Foundation GACR 20-21727X. The
code has been originally developed by L. Papa, and subsequently optimized by
M. Corrias.

## References

  * Giessibl (2003) F. J. Giessibl, Rev. Mod. Phys. 75, 949 (2003).
  * Morita  _et al._ (2015) S. Morita, F. Giessibl, E. Meyer, and R. Wiesendanger,  _Noncontact Atomic Force Microscopy: Volume 3_, NanoScience and Technology (Springer International Publishing, 2015).
  * Pavliček and Gross (2017) N. Pavliček and L. Gross, Nature Reviews Chemistry 1, 1 (2017).
  * Scheuerer  _et al._ (2019) P. Scheuerer, L. L. Patera, F. Simbürger, F. Queck, I. Swart, B. Schuler, L. Gross, N. Moll, and J. Repp, Phys. Rev. Lett. 123, 066001 (2019).
  * Kempkes  _et al._ (2019) S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. M. Smith, Nature physics 15, 127 (2019).
  * Fournier  _et al._ (2011) N. Fournier, C. Wagner, C. Weiss, R. Temirov, and F. S. Tautz, Phys. Rev. B 84, 035435 (2011).
  * Wagner  _et al._ (2012) C. Wagner, N. Fournier, F. S. Tautz, and R. Temirov, Phys. Rev. Lett. 109, 076102 (2012).
  * Gross  _et al._ (2009) L. Gross, F. Mohn, N. Moll, P. Liljeroth, and G. Meyer, Science 325, 1110 (2009).
  * Alldritt  _et al._ (2020) B. Alldritt, P. Hapala, N. Oinonen, F. Urtev, O. Krejci, F. F. Canova, J. Kannala, F. Schulz, P. Liljeroth, and A. S. Foster, Science Advances 6, eaay6913 (2020).
  * Krull  _et al._ (2020) A. Krull, P. Hirsch, C. Rother, A. Schiffrin, and C. Krull, Communications Physics 3, 1 (2020).
  * Kalinin  _et al._ (2021a) S. V. Kalinin, M. Ziatdinov, J. Hinkle, S. Jesse, A. Ghosh, K. P. Kelley, A. R. Lupini, B. G. Sumpter, and R. K. Vasudevan, ACS Nano 15, 12604 (2021a).
  * Kalinin  _et al._ (2022) S. V. Kalinin, C. Ophus, P. M. Voyles, R. Erni, D. Kepaptsoglou, V. Grillo, A. R. Lupini, M. P. Oxley, E. Schwenker, M. K. Y. Chan, J. Etheridge, X. Li, G. G. D. Han, M. Ziatdinov, N. Shibata, and S. J. Pennycook, Nature Reviews Methods Primers 2, 11 (2022).
  * Choudhary  _et al._ (2022) K. Choudhary, B. DeCost, C. Chen, A. Jain, F. Tavazza, R. Cohn, C. W. Park, A. Choudhary, A. Agrawal, S. J. L. Billinge, E. Holm, S. P. Ong, and C. Wolverton, npj Computational Materials 8, 59 (2022).
  * Rashidi and Wolkow (2018) M. Rashidi and R. Wolkow, ACS Nano 12 (2018).
  * Gordon  _et al._ (2020) O. M. Gordon, F. L. Junqueira, and P. J. Moriarty, Machine Learning: Science and Technology 1, 015001 (2020).
  * Ziatdinov  _et al._ (2018) M. Ziatdinov, O. Dyck, A. Maksov, B. M. Hudak, A. R. Lupini, J. Song, P. C. Snijders, R. K. Vasudevan, S. Jesse, and S. V. Kalinin, arXiv e-prints (2018), arXiv:1801.05133 .
  * Kalinin  _et al._ (2021b) S. V. Kalinin, M. P. Oxley, and M. Valleti, npj Comput Mater 7, 181 (2021b).
  * Ziletti  _et al._ (2018) A. Ziletti, D. Kumar, M. Scheffler, and L. M. Ghiringhelli, Nature Communications 9, 2775 (2018).
  * Leitherer  _et al._ (2021) A. Leitherer, A. Ziletti, and L. M. Ghiringhelli, Nature Communications 12, 6234 (2021).
  * Carracedo-Cosme  _et al._ (2021) J. Carracedo-Cosme, C. Romero-Muñiz, and R. Pérez, Nanomaterials 11 (2021), 10.3390/nano11071658.
  * Kalinin  _et al._ (2016) S. V. Kalinin, E. Strelcov, A. Belianinov, S. Somnath, R. K. Vasudevan, E. J. Lingerfelt, R. K. Archibald, C. Chen, R. Proksch, N. Laanait, and S. Jesse, ACS Nano 10, 9068 (2016).
  * Aversa  _et al._ (2018) R. Aversa, M. H. Modarres, S. Cozzini, R. Ciancio, and A. Chiusole, Scientific Data 5, 180172 (2018).
  * Schwenker  _et al._ (2020) E. Schwenker, F. Sen, C. Wolverton, C. Ophus, and M. K. Chan,  (2020), 10.18126/SZEQ-YDE5.
  * Ede (2020) J. M. Ede, Machine Learning: Science and Technology 1, 045003 (2020).
  * Choudhary  _et al._ (2019) K. Choudhary, K. F. Garrity, C. Camp, S. V. Kalinin, R. Vasudevan, M. Ziatdinov, and F. Tavazza,  (2019), 10.48550/ARXIV.1912.09027.
  * Carracedo-Cosme  _et al._ (2022) J. Carracedo-Cosme, C. Romero-Muñiz, P. Pou, and R. Pérez, Journal of Chemical Information and Modeling 62, 1214 (2022), pMID: 35234034, https://doi.org/10.1021/acs.jcim.1c01323 .
  * Somnath  _et al._ (2019) S. Somnath, C. R. Smith, N. Laanait, R. K. Vasudevan, and S. Jesse, Microscopy and Microanalysis 25, 220–221 (2019).
  * Ziatdinov and Kalinin (2021) M. Ziatdinov and S. Kalinin, Microscopy and Microanalysis 27, 3000–3002 (2021).
  * Nord  _et al._ (2017) M. Nord, P. Vullum, and I. MacLaren, Adv Struct Chem Imag 3, 9 (2017).
  * Belianinov  _et al._ (2015) A. Belianinov, Q. He, M. Kravchenko, S. Jesse, A. Borisevich, and S. V. Kalinin, Nature Communications 6, 7801 (2015).
  * Somnath  _et al._ (2018) S. Somnath, C. R. Smith, S. V. Kalinin, M. Chi, A. Borisevich, N. Cross, G. Duscher, and S. Jesse, Advanced Structural and Chemical Imaging (2018).
  * Laanait  _et al._ (2016) N. Laanait, M. Ziatdinov, Q. He, and A. Borisevich, Advanced Structural and Chemical Imaging 2, 14 (2016).
  * Lowe (2004) D. G. Lowe, International Journal of Computer Vision 60, 91 (2004).
  * Gan  _et al._ (2007) G. Gan, C. Ma, and J. Wu,  _Data clustering: theory, algorithms, and applications_, Vol. 20 (Siam, 2007).
  * Hapala  _et al._ (2014) P. Hapala, G. Kichin, C. Wagner, F. S. Tautz, R. Temirov, and P. Jelinek, Physical Review B 90, 085421 (2014).
  * Hapala  _et al._ (2016) P. Hapala, M. Švec, O. Stetsovych, N. J. van der Heijden, M. Ondráček, J. van der Lit, P. Mutombo, I. Swart, and P. Jelínek, Nature Communications 7, 1 (2016).
  * Bradski (2000) G. Bradski, Dr. Dobb’s Journal of Software Tools (2000).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).
  * Rousseeuw (1987) P. J. Rousseeuw, Journal of Computational and Applied Mathematics 20, 53 (1987).
  * Hommrich  _et al._ (2002) J. Hommrich, S. Hümann, and K. Wandelt, Faraday Discuss. 121, 129 (2002).
  * Lee  _et al._ (2006) S. Lee, J. Park, R. Ragan, S. Kim, Z. Lee, D. K. Lim, D. A. A. Ohlberg, and R. S. Williams, Journal of the American Chemical Society 128, 5745 (2006).
  * Bai and Breen (2008) L. Bai and D. Breen, Journal of Graphics Tools 13, 53 (2008).
  * Giessibl (2019) F. J. Giessibl, Review of Scientific Instruments 90, 011101 (2019).
  * Sokolović  _et al._ (2019) I. Sokolović, M. Schmid, U. Diebold, and M. Setvin, Phys. Rev. Materials 3, 034407 (2019).
  * Sokolović  _et al._ (2021) I. Sokolović, G. Franceschi, Z. Wang, J. Xu, J. Pavelec, M. Riva, M. Schmid, U. Diebold, and M. Setvín, Physical Review B 103, L241406 (2021).
  * Sokolović  _et al._ (2020) I. Sokolović, M. Reticcioli, M. Čalkovskỳ, M. Wagner, M. Schmid, C. Franchini, U. Diebold, and M. Setvín, Proceedings of the National Academy of Sciences 117, 14827 (2020).
  * Reticcioli  _et al._ (2017) M. Reticcioli, M. Setvín, X. Hao, P. Flauger, G. Kresse, M. Schmid, U. Diebold, and C. Franchini, Physical Review X 7, 031053 (2017).
  * Reticcioli  _et al._ (2019) M. Reticcioli, I. Sokolović, M. Schmid, U. Diebold, M. Setvin, and C. Franchini, Phys. Rev. Lett. 122, 016805 (2019).
  * Meier  _et al._ (2022) M. Meier, J. Hulva, Z. Jakub, F. Kraushofer, M. Bobić, R. Bliem, M. Setvin, M. Schmid, U. Diebold, C. Franchini, and G. S. Parkinson, Science Advances 8, eabn4580 (2022).
  * Kunat  _et al._ (2009) M. Kunat, F. Traeger, D. Silber, H. Qiu, Y. Wang, A. C. van Veen, C. Wöll, P. M. Kowalski, B. Meyer, C. Hättig, and D. Marx, The Journal of Chemical Physics 130, 144703 (2009).
  * Prates Ramalho  _et al._ (2017) J. P. Prates Ramalho, F. Illas, and J. R. B. Gomes, Phys. Chem. Chem. Phys. 19, 2487 (2017).