# YASTN: Yet another symmetric tensor networks;
A Python library for abelian symmetric tensor network calculations.

Marek M. Rams  Jagiellonian University, Institute of Theoretical Physics,
Łojasiewicza 11, 30-348 Kraków, Poland<EMAIL_ADDRESS>Gabriela
Wójtowicz  Jagiellonian University, Institute of Theoretical Physics,
Łojasiewicza 11, 30-348 Kraków, Poland  Institut für Theoretische Physik und
IQST, Albert-Einstein-Allee 11, Universität Ulm, D-89081 Ulm, Germany
<EMAIL_ADDRESS>Aritra Sinha  Jagiellonian University, Institute
of Theoretical Physics, Łojasiewicza 11, 30-348 Kraków, Poland  Max Planck
Institute for the Physics of Complex Systems, Nöthnitzer Strasse 38, Dresden
01187, Germany<EMAIL_ADDRESS>Juraj Hasik  Institute for Theoretical
Physics and Delta Institute for Theoretical Physics, University of Amsterdam,
Science Park 904, 1098 XH Amsterdam, The Netherlands  Department of Physics,
University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
<EMAIL_ADDRESS>
###### Abstract

We present an open-source tensor network Python library for quantum many-body
simulations. At its core is an abelian-symmetric tensor, implemented as a
sparse block structure managed by logical layer on top of dense multi-
dimensional array backend. This serves as the basis for higher-level tensor
networks algorithms, operating on matrix product states and projected
entangled pair states, implemented here. Using appropriate backend, such as
PyTorch, gives direct access to automatic differentiation (AD) for cost-
function gradient calculations and execution on GPUs or other supported
accelerators. We show the library performance in simulations with infinite
projected entangled-pair states, such as finding the ground states with AD, or
simulating thermal states of the Hubbard model via imaginary time evolution.
We quantify sources of performance gains in those challenging examples allowed
by utilizing symmetries.

##  I Introduction

Full numerical treatment of quantum-mechanical systems is generally
prohibitively expensive due to the exponential growth of Hilbert space size
with the number of interacting degrees of freedom. Tensor network (TN)
techniques allow for the efficient representation and manipulation of states
of such large quantum systems [1, 2, 3]. The density matrix renormalization
group (DMRG) introduced by White [4, 5] and its modern reformulation in terms
of matrix product states [6, 7, 8, 9, 10] (MPS), a one-dimensional tensor
network, is a prime example of TN capabilities. Since their inception, MPS
quickly became a reference method for addressing ground states in one
dimension, and were soon followed by extensions to excited states, time
evolution, and open systems forming a comprehensive framework.

The descriptive power of tensor networks is not limited to systems in one
dimension. Despite their intrinsic 1D geometry, MPS can be readily applied to
models in higher dimensions by linearly ordering the sites. Typically, for 2D
models a finite cylinder is chosen and mapped to the MPS ansatz by imposing
ordering winding around the cylinder circumference as shown in Fig. 1(c). More
natural TN geometry for two-dimensional states is assured by the projected
entangled-pair states [11, 12] (PEPS) (see Fig. 1e) and the similar for 3D
states [13, 14]. The TN ansätze in Fig. 1, provide a state-of-the-art
numerical approach to strongly correlated systems of condensed matter. The
computational complexity of MPS typically scales as $\mathcal{O}(D^{3})$, and
PEPS algorithms often reach $\mathcal{O}(D^{12})$ scaling, where bond
dimension $D$ governs the size of the tensors and the overall precision of the
TN approximations. While the scaling of PEPS seems less favorable, it is
important to note that the bond dimension encodes correlations between sites.
Under the snake-MPS ansatz, the correlations inside a column are stretched due
to the row-by-row mapping, leading PEPS to reach comparable or better
precision even at low bond dimensions once the cylinders become too wide.

Figure 1: Tensor networks. Diagrams depict (a) a rank-3 tensor of the MPS with
left and right virtual spaces $V$ and physical space $\mathcal{H}$, (b) an MPS
ansatz, (c) a MPS winding on a finite width-6 cylinder, (d) a rank-5 tensor of
the PEPS, and (e) an infinite-PEPS ansatz.

The most effective way to mitigate the computational complexity is to leverage
the symmetries present in physical systems. Two principal types of symmetries
to consider are spatial and internal symmetries. Tensor networks can be
formulated directly in the thermodynamic limit by infinitely repeated pattern
(a unit cell) of tensors, hence realizing translation symmetry. These are
infinite-MPS (iMPS) also known as uniform MPS [15] in one dimension and the
infinite-PEPS (iPEPS) [16] in two dimensions. The computational complexity of
iMPS/iPEPS algorithms scales linearly with the size of the unit cell. For
internal symmetries $U|\Psi\rangle=|\Psi\rangle$, we consider their common
form of global symmetries, i.e., when $U=\otimes_{i}u_{i}$ with the same
unitaries $u_{i}$ acting on each site. These can be both abelian (e.g.,
particle conservation) or non-abelian (e.g., SU(2)-spin). Crucially, such
global symmetries can be implemented in TNs locally, by requiring individual
tensors to transform covariantly under the action $u$ of the symmetry group
[17, 18, 19, 20, 21]. These symmetric tensors take block-spare form, with
original dense virtual spaces $V$ of bond dimension $D$ split into a direct
sum of blocks $V=\oplus_{r}V_{r}$ with dimensions $\\{D_{1},\ldots,D_{r}\\}$
each associated to irreducible representation $r$ of the considered symmetry
group. The use of block sparsity substantially lowers computational complexity
permitting large-$D$ simulations, in particular for (i)PEPS algorithms.

Here, we introduce the Yet Another Symmetric Tensor Network (YASTN) library
[22]. YASTN is an open-source Python library with abelian-symmetric tensor as
a basic type and associated linear algebra operations on such tensors. The
implementation enables automatic differentiation (AD) via appropriate dense
linear algebra backends, allowing convenient variational optimization of TNs.
This is particularly important for iPEPS [23, 24, 25], where no alternative
direct energy minimization algorithms are known. This is in contrast to the
(i)MPS where the DMRG provides efficient and robust optimization. YASTN thus
joins a continually growing collection of tensor network software with various
degree of support for symmetries and automatic differentiation such as iTensor
[26], TenPy [27], Block2 [28], Quantum TEA [29], TensorNetwork [30], Cytnx
[31], TeNes [32], TensorKit [33], Qspace [34], peps-torch [35], ad-peps [36],
variPEPS [58], PEPSKit [59].

In the following sections, we outline the design principles of YASTN and
present a set of benchmarks demonstrating the computational speed-up from
abelian symmetries. We focus on variational optimization of iPEPS for
SU(2)-symmetric spin-$\frac{1}{2}$ model, SU(3)-symmetric model, and
observables of a Hubbard model at finite temperature simulated via imaginary-
time evolution.

##  II Design principles

In this section, we give an overview of the structure of YASTN, presented in
Fig. 2, and comment on some aspects of implementation. The basic building
block of the library is the yastn.Tensor which is defined by the symmetry
structure and the backend. The symmetry structure determines allowed blocks
and how to manipulate them when performing tensor algebra. The backend handles
the execution of dense linear algebra operations and storage of tensor
elements. These two are independent of each other. Symmetric tensors are used
to construct TN ansätze such as MPS and PEPS, and finally define high-level
algorithms that are applied to specific TN. For a detailed description of the
library and all its functionalities, see the documentation under Ref. [22].

Figure 2: Schematic design of yastn package. The core element of the package
is yastn.Tensor, which implements block-sparse tensor structure corresponding
to a given abelian symmetry on top of a dense linear algebra backend, such as
NumPy [38] or PyTorch [39]. More complex tensor networks, built on symmetric
tensor, include standard but versatile MPS toolbox and 2D PEPS implementations
to simulate the time evolution or ground-state variational optimization. The
latter can benefit from automatic differentiation supported by some underlying
backends, such as PyTorch.

###  II.1 Abelian-symmetric tensor

Tensors are multilinear maps from products of several vector spaces

$T:\quad V^{1}\otimes V^{2}\otimes V^{3}\otimes...\rightarrow\mathbb{C},$ (1)

where $V^{i}$ is a vector space and $\otimes$ is a tensor product. In quantum
mechanical context we work with Hilbert spaces $\mathcal{H}$ and their duals
$\mathcal{H}^{*}$, and by choosing some bases in each of these spaces the
tensors can be written out in components

$T=\sum_{abc...ijk...}T^{abc...}_{ijk...}|i\rangle|j\rangle|k\rangle\ldots\langle
a|\langle b|\langle c|\ldots,$ (2)

where $i,j,k\dots$ are indices of bases in $\mathcal{H}$ spaces, $a,b,c\dots$
in $\mathcal{H}^{*}$ spaces, and $T^{abc...}_{ijk...}$ is the corresponding
tensor element. The action of the element $g$ of abelian group $G$ on tensor
$T$ can be represented in a proper basis by diagonal matrices $U^{i}(g)$
transforming the tensor elements

$\begin{split}(g\circ
T)^{ab...}_{ij...}=&\sum_{a^{\prime}b^{\prime}...i^{\prime}j^{\prime}...}T^{a^{\prime}b^{\prime}...}_{i^{\prime}j^{\prime}...}{[U^{1}(g)}]^{i^{\prime}}_{i}{[U^{2}(g)}]^{j^{\prime}}_{j}\ldots\\\
&\quad\times[U^{m}(g)^{*}]^{a}_{a^{\prime}}[U^{m+1}(g)^{*}]^{b}_{b^{\prime}}\ldots.\end{split}$
(3)

The matrix elements of $U^{i}(g)$ are

${[U^{i}(g)]}^{j}_{j^{\prime}}=\delta_{jj^{\prime}}\,\textrm{exp}({-i\theta_{g}t^{[i]}_{j}}),$
(4)

forming a diagonal matrix of complex phases defined by integer charges
$t^{[i]}_{j}$, with angle $\theta_{g}\in[0,2\pi)$ which depends on $g\in G$,
and $\delta_{jj^{\prime}}$ being Kronecker delta. Therefore, under the action
of $g\in G$ each tensor element simply acquires a phase given by the sum of
charges

$\begin{split}(g\circ
T)^{ab...}_{ij...}=T^{ab...}_{ij...}\textrm{exp}[-i\theta_{g}(&t^{[1]}_{i}+t^{[2]}_{j}+\ldots\\\
&-t^{[m]}_{a}-t^{[m+1]}_{b}-\ldots)].\end{split}$ (5)

This form of the transformation gives a simple selection rule, a charge
conservation, on the elements of symmetric tensors

$t^{[1]}_{i}+t^{[2]}_{j}+...-t^{[m]}_{a}-t^{[m+1]}_{b}-...=n.$ (6)

The charge of each non-zero element $T^{ab...}_{ij...}$ of a symmetric tensor
must be $n$. In the case of $n=0$, such tensor is invariant (unchanged) under
the action of the symmetry. Otherwise, it transforms covariantly as all its
elements are altered by the same complex phase $\textrm{exp}(-i\theta_{g}n)$.
The charges $t^{[i]}_{j}$ and $n$ and the precise form of their addition
depend on the abelian group. For elementary abelian groups such as $Z_{n}$ or
$U(1)$ the individual charges $t^{[i]}_{j}$ are elements of $\mathbb{Z}_{n}$
or $\mathbb{Z}$ respectively, while for direct products of abelian groups,
they become vectors in corresponding product of $\mathbb{Z}_{n}$’s and
$\mathbb{Z}$’s.

By ordering the basis elements in each Hilbert space by their charge, the
tensor $T^{ab...}_{ij...}$ naturally attains a block-sparse structure which is
central to the computational advantage offered by abelian-symmetric tensor
network algorithms.

At the core of YASTN is the implementation of symmetric tensor yastn.Tensor,
as outlined in Refs. [18, 19, 17]. It is defined jointly by symmetry (block)
structure data and tensor elements of existing blocks. First, we define a
vector space with a charge structure, a yastn.Leg, determined by a signature
$s=\pm 1$ (distinguishing between $\mathcal{H}$ and dual $\mathcal{H}^{*}$),
its charge sectors $\bf t$, and their corresponding dimensions $\bf D$, now in
boldface to underline their vector nature,

$V(s,{\bf t},{\bf D})=\oplus_{\rho\in{\bf t}}\mathbb{C}^{{\bf D}_{\rho}},$ (7)

where $\rho$ now enumerates different charges instead of basis 111The
conjugation of the leg, i.e., mapping space $\mathcal{H}$ to its dual space
$\mathcal{H}^{*}$, is equivalent to flip of the signature and complex
conjugation of elements.. This space is a direct sum of simple spaces
$\mathbb{C}^{{\bf D}_{\rho}}$, dubbed charge sectors. The effective dimension
of such space is the sum of dimensions of individual charge sectors

$D=\sum_{\rho}{\bf D}_{\rho}.$ (8)

In the remainder of the text, we will refer to $D$ as the bond dimension, when
discussing scaling of computational complexity or memory requirements of TN
algorithms with symmetric tensors.

The abelian symmetric tensor of rank-$N$ is specified by the product of $N$
such vector spaces

$T:\otimes^{N}_{i=1}V(s^{[i]},{\bf t}^{[i]},{\bf
D}^{[i]})\rightarrow\mathbb{C},$ (9)

The following is an example creating a random $U(1)$-symmetric tensor with
total charge 222One can always define tensors with an extra dummy leg
$V(-1,(n,),(1,))$, having a single charge sector of a unit dimension, making
it invariant under symmetry transformations. $n=0$ with specified legs:

⬇

1import yastn

2from yastn.backend import backend_np

3from yastn.sym import sym_U1

4

5u1 = yastn.make_config(sym=sym_U1, backend=backend_np)

6l = yastn.Leg(u1, s=1, t=(-1,1), D=(1,1))

7lc = l.conj()

8H = yastn.rand(u1, legs=[l,l,lc,lc], n=0)

which is, for example, compatible with a Hamiltonian
$H=\vec{S}_{1}\cdot\vec{S}_{2}$ of two spin-$\frac{1}{2}$ degrees of freedom.
The configuration created by yastn.make_config specifies the symmetry, e.g.,
yastn.sym.sym_U1 for the $U(1)$ in the example, and dense linear algebra
backend (see below), e.g., yastn.backend.backend_np for NumPy. The covariant
transformation property of $T$ is imposed by the charge conservation of non-
zero blocks. Any block of tensor $T$ can be identified by selecting a charge
sector $\rho_{i}\in{\bf t}^{[i]}$ on each of the legs, i.e., an $N$-tuple of
charges $(\rho_{0},\rho_{1},...,\rho_{N})$. All non-zero blocks must satisfy

$\sum_{i=1}^{N}s^{[i]}\rho_{i}=n,$ (10)

which is the block-sparse version of element-wise charge conservation rule of
Eq. 6.

Finally, we remark on the storage of tensor elements. In YASTN, the block data
is initialized lazily. The storage is allocated only for the blocks which have
been assigned a non-zero value, i.e., blocks allowed by the charge
conservation but not assigned any value are not stored. All allocated blocks
are serialized together in a 1D array.

###  II.2 Fusion and contractions

The key operations on symmetric tensors, necessary for manipulating tensor
networks, are tensor reshape and permutation, commonly dubbed fusion in this
context, and tensor contractions. Fusion resolves the tensor product of
several spaces as a new space, i.e., fusion of legs into a new leg. Unlike
reshaping of the dense tensor, the shape cannot be chosen freely. Instead, it
is determined by the structure of fused spaces. In particular, fusion orders
and accumulates tensor products of charge sectors on selected legs into new
charge sectors under the joint leg

$V(s^{[i]},{\bf t}^{[i]},{\bf D}^{[i]})\otimes V(s^{[j]},{\bf t}^{[j]},{\bf
D}^{[j]})\rightarrow V(s^{[r]},{\bf t}^{[r]},{\bf D}^{[r]}),$ (11)

with new charge sectors ${\bf t}^{[r]}$ given by the unique combinations of
charges ${\bf t}^{[i]}\otimes{\bf t}^{[j]}$

${\bf t}^{[r]}:=\\{\nu=s^{[r]}(s^{[i]}\rho+s^{[j]}\rho^{\prime}):\rho\in{\bf
t}^{[i]},\rho^{\prime}\in{\bf t}^{[j]}\\}.$ (12)

The dimension of new charge sector $\nu\in{\bf t}^{[r]}$ is 333Following a
lazy approach adopted in YASTN, a new fused leg contains only charges for
which some non-zero tensor block exists. As such, fusion in YASTN is always
done in the context of particular tensor.

${\bf D}^{[r]}_{\nu}=\sum_{\begin{subarray}{c}\rho,\rho^{\prime}\\\
\nu=s^{[r]}(s^{[i]}\rho+s^{[j]}\rho^{\prime})\end{subarray}}{\bf
D}^{[i]}_{\rho}{\bf D}^{[j]}_{\rho^{\prime}}.$ (13)

The fusion and un-fusion calls are demonstrated below on previously
constructed rank-4 tensor $H$, first fusing pairs of legs resulting in a
matrix form

⬇

1H_mat = H.fuse_legs(axes=((0,1), (2,3)))

2H = H_mat.unfuse_legs(axes=(0,1))

In the example above, the YASTN first computes the structure of the resulting
tensor with fused leg, i.e., the tuple $(s^{[r]},{\bf t}^{[r]},{\bf D}^{[r]})$
and a set of dense linear algebra jobs (permutes, reshapes, and copies) to be
executed by the backend to populate the new 1D storage array with serialized
blocks. The resulting tensor records the original structure and hence, the
fusion can be reverted to restore the original structure.

The tensor contraction of symmetric tensors is realized by a commonly adopted
workflow. First, the tensors are fused into matrices 444For valid contraction,
the structures of the contracted legs must be compatible, including the
origins of any fused leg. YASTN automatically resolves a situation when some
charge in the fusion history is missing in a to-be-contracted fused leg but is
present in its contraction partner, utilizing information on the tensor’s
fusion history stored in each tensor., then multiplied along the contracted
legs, and finally unfused to obtain the desired form:

$\begin{split}\sum_{\\{x\\}}A_{\\{i\\}\cup\\{x\\}}B_{\\{j\\}\cup\\{x\\}}&\overset{\text{fuse}}{\longrightarrow}\sum_{X}A_{IX}B_{XJ}=\\\
C_{IJ}&\overset{\text{unfuse}}{\longrightarrow}C_{\\{i\\}\\{j\\}}\end{split}$
(14)

where $\\{x\\}$ is a set of common legs that become fused into single leg $X$,
and original legs $\\{i\\}$ and $\\{j\\}$ are restored from the fused $I$ and
$J$ to obtain the final tensor. Here, we first show an example call for
pairwise tensor contraction, and second, an equivalent given in terms of
explicit operations

⬇

1H2 = yastn.tensordot(H, H, axes=((2,3), (0,1)))

2H2 = (H_mat @ H_mat).unfuse_legs(axes=(0,1))

For both fusion and multiplication, the YASTN first precomputes what are the
non–zero blocks so the backend performs only the relevant operations. Using
these elementary operations, contractions of more general networks are
supported through convenience functions, such as the einsum and ncon [40]
functions (in our case differing only by syntax).

###  II.3 Tensor network algorithms

The symmetric tensor serves as a foundation for higher-level tensor network
structures and algorithms. Here, the YASTN comes with MPS and PEPS modules.
The MPS module supports finite-size MPS with the implementation of a range of
standard algorithms, including DMRG for ground-state optimization, TDVP [41]
for time evolution, and the overlap maximization [42] against a general
target, i.e., MPS, sum of MPSs, or sum of MPO-MPS products. This is
complemented by a versatile high-level (Hamiltonian) MPO generator. The MPS
module provides subroutines for some PEPS methods, e.g., it was utilized in
[43] for boundary MPS contraction and long-range correlations calculation in a
finite PEPS defined on a cylinder. At the same time, it is a versatile
computational toolbox on its own. For instance, it has been employed in
simulations of Lindbladian dynamics in the context of fermionic quantum
transport [44], where the $U(1)$ symmetry reflects a lack of correlations
between different particle-number sectors of a density matrix.

The PEPS module features the implementation of fermionic PEPS, dubbed fPEPS
(which also allows simulations of systems without fermionic statistics). It
covers both finite PEPS defined on a square grid and its infinite versions for
translationally invariant (over a unit-cell) systems in the thermodynamic
limit. It supports a range of time-evolution algorithms, starting with
neighborhood tensor update (NTU) scheme [45, 46], its refinement to a family
of larger environmental clusters [43], ending on a full-update type of schemes
[16, 47]. It is viable for imaginary-time evolution, e.g., in the context of
finite temperature simulation of density matrix purification [48] or
minimally-entangled typical thermal states [49], and real-time simulations,
e.g., of pure state quench-dynamics in disordered spin systems [43].

##  III Examples

To demonstrate the use and the versatility of YASTN we present three end-to-
end numerical examples centered on iPEPS. We show computational speed-up and
reduced memory footprint obtained with YASTN by utilizing abelian symmetries
for the following examples:

  1. 1.

Sec. III.1: variational optimization of iPEPS with $D\leq 8$ for
antiferromagnetic spin-$\frac{1}{2}$ model on a square lattice using $U(1)$
symmetry,

  2. 2.

Sec. III.2: variational optimization of iPEPS with $D\leq 13$ for SU(3) model
on Kagome lattice using $U(1){\times}U(1)$ symmetry,

  3. 3.

Sec. III.3: observables of thermal iPEPS of Hubbard model at finite
temperature using $Z_{2}$, $U(1)$, and $U(1){\times}U(1)$ symmetry, with $D$
up to $36$ for the latter.

In Sec. III.1 and Sec. III.2 we optimize iPEPS for SU(2) model on square and
SU(3) model on Kagome lattices respectively. First, given an iPEPS generated
by a set of tensors $\vec{a}=\\{a,b,\ldots\\}$, we compute an approximate
environment tensors $\vec{E}(\vec{a})$ (specified below) with the precision
governed by the environment dimension $\chi$. Then, environment $\vec{E}$ and
tensors $\vec{a}$ are combined to evaluate the energy per site $e$ of the
Hamiltonian. Finally, the reverse mode of AD (i.e., backpropagation) is
invoked to compute the gradient $\partial e/\partial\vec{a}$. The most
computationally intensive stage is the construction of the environments,
scaling as the cube of $D^{2}\chi$, which assuming the necessary $\chi\propto
D^{2}$ gives the overall complexity $\mathcal{O}(D^{12})$, where $D$ is the
bond dimension of iPEPS tensors. We use iPEPS optimization implemented in
peps-torch [35], here operating on YASTN’s symmetric tensors.

Figure 3: Corner transfer matrix iteration. In (a), CTM approximates parts of
an infinite tensor network by a set of finite environment tensors
characterized by environment bond dimension $\chi$. We depict elements of the
CTM algorithm step (iteration) that dominate the computational effort shown in
Figs.4–6. Panels (b), (b’), and (c) depict the construction of an enlarged
corner combining CTM environment tensors (rectangular) with PEPS tensors
(circle). Panel (d) shows an SVD decomposition of a product of four enlarged
corners that is then used to construct the CTM projections from enlarged
virtual spaces.

For the presented examples, we employ the corner transfer matrix (CTM)
algorithm [50, 51, 52] to compute the environments. CTM approximates
environments $\vec{E}=\\{C,T\\}$ of iPEPS by a set of $\chi\times\chi$ corner
matrices $C$ and $\chi\times D^{2}\times\chi$ transfer tensors $T$, as shown
in Fig. 3(a). Alternatively, one can use boundary MPS methods [16, 42, 53].
The computational complexity of CTM arises from two sources, see Fig. 3,
tensor contractions and singular value decomposition (SVD) when computing low-
rank approximations, both scaling as $\mathcal{O}(D^{12})$. In practice, for
simulations without symmetries the SVD gives a substantially greater
contribution due to a higher scaling prefactor and a poor speed-up offered by
multithreading or GPU acceleration compared to tensor contractions. However,
for symmetric iPEPS the situation becomes more nuanced as we demonstrate in
the examples below.

In Sec. III.3, the purification techniques are applied to compute thermal
expectation values for the Fermi-Hubbard model on a $2$D square lattice. The
ancilla trick is used to effectively transform the thermal density matrix into
a purified wavefunction and then imaginary time evolution is performed to
reach the target temperature. We adopt the NTU algorithm to optimize the time-
evolution, with computational scaling of $\mathcal{O}(D^{8})$ dominated by
tensor contractions. We focus our example on the final calculation of the
expectation values using CTM with $\chi=5D$, translating to
$\mathcal{O}(D^{9})$ scaling. We quantify the sources of advantage offered by
incorporation of higher symmetries.

###  III.1 Heisenberg antiferromagnet with anisotropy

We revisit a Heisenberg model with anisotropy in the couplings describing a
system of coupled spin-$\frac{1}{2}$ ladders

${\cal H}=J\sum_{R}{\bf S}_{R}\cdot{\bf S}_{R+\hat{x}}+\sum_{R}J_{R}{\bf
S}_{R}\cdot{\bf S}_{R+\hat{y}},$ (15)

where ${\bf S}_{R}=(S^{x}_{R},S^{y}_{R},S^{z}_{R})$ is the $S=\frac{1}{2}$
operator on the site $R=(x,y)$ of a square lattice, $\hat{x}$ and ${\hat{y}}$
are unit vectors in $x$ and $y$ direction, and $J_{R}=J$ or $J_{R}=\alpha J$,
depending on the parity of $y$. By varying $\alpha$, this model interpolates
between the Heisenberg model on the square lattice at $\alpha=1$ and a system
of decoupled two-leg ladders at $\alpha=0$. The model was previously addressed
by iPEPS in Ref. [54]. We adopt the same description that uses four non-
equivalent tensors $\vec{a}=\\{a,b,c,d\\}$ arranged in a $2{\times}2$ unit
cell 555A more efficient description might generate all tensors in
$2{\times}2$ unit cell from a single parent tensor $a$ by use of unitary
$-i\sigma^{y}$ acting on physical index and/or permutation of virtual indices
generated by the reflection along the $x$-axis. Nevertheless, such
parametrization would not change CTM complexity..

The ground states possess $U(1)$-symmetry corresponding to the conservation of
$S^{z}$ which we leverage to work with $U(1)$-symmetric iPEPS. The results,
summarized in Fig. 4, show a rapidly growing computational advantage of
symmetric iPEPS for bond dimensions $D>4$. While at $D=4$ the overhead due to
the block-sparse logic is still significant, at the largest bond dimension
considered, $D=8$, a $30$-fold speed-up is observed. In practical terms, the
convergence of the CTM towards the desired precision, here measured by the
error on the energy per site becoming lower than $\epsilon<10^{-8}$, typically
requires $\mathcal{O}(10)$ iterations, thus without $U(1)$ symmetry a single
optimization step would already take hours. Details of the block-sparse
structure and its impact on the CTM are visualized in Fig. 4(b,c). At the
largest bond dimension, $D=8$, the fusion of enlarged corner into a block-
diagonal matrix requires processing of roughly $\mathcal{O}(100)$ blocks by
performing dense permutes, reshapes, and copies, with the largest block having
$\mathcal{O}(10^{5})$ elements. The cost of subsequent SVD is dominated by the
largest block(s) of fused enlarged corner, which are $L\times L$ matrices with
$L=1500\sim 2000$. As a result, the computational time contributions are
shared between the SVD and tensor contractions with fusions roughly as 3:2,
with SVD being the dominant factor.

Figure 4: Use case 1: optimization of $U(1)$-symmetric iPEPS for model of
spin-$\frac{1}{2}$ coupled ladders. In (a), we show the scaling of the wall
time per CTM step (in seconds) for the entire gradient optimization step of
iPEPS. The bond dimension of the environment is $\chi=2D^{2}$. The inset shows
the relative speed-up compared to an implementation without symmetry. In (b),
we show a distribution of blocks of an enlarged corner by their size (number
of elements) before fusion to a block-diagonal matrix as shown in Fig. 3(c).
In (c), we show the sizes of $L\times L$ blocks after the fusion.

###  III.2 SU(3) model on Kagome lattice

We consider an SU(3)-symmetric model on Kagome lattice, analyzed recently in
Ref. [55], where each site holds a single degree of freedom from the
fundamental representation $\mathbf{3}$ of SU(3) group spanned by three states
$\\{|\alpha\rangle,|\beta\rangle,|\gamma\rangle\\}$. The Hamiltonian reads

$\displaystyle H$ $\displaystyle=$ $\displaystyle J\sum_{\langle
i,j\rangle}P_{ij}+\sum_{\triangle ijk}(KP_{ijk}+h.c.),$ (16)

where $P_{ij}$ is a permutation of local states on nearest-neighbour bonds
such that
$P_{ij}|\alpha\rangle_{i}|\beta\rangle_{j}=|\beta\rangle_{i}|\alpha\rangle_{j}$,
$P_{ijk}$ is a clockwise permutation of local states on nearest-neighbour
triangles such that
$P_{ijk}|\alpha\rangle_{i}|\beta\rangle_{j}|\gamma\rangle_{k}=|\gamma\rangle_{i}|\alpha\rangle_{j}|\beta\rangle_{k}$,
with fixed choice of the orientation of triangles, and $J$ and $K$ are real-
and complex-valued couplings respectively.

In this section, we demonstrate an advantage of $U(1){\times}U(1)$-symmetric
iPEPS, utilizing maximal abelian subgroup of SU(3), over implementation
without symmetries 666In this example, besides iPEPS, one can also use
different ways to construct two-dimensional TN ansatz on Kagome lattice, i.e.,
infinite projected simplex states (iPESS), however, the computational
complexity $\mathcal{O}(D^{12})$, attributable to CTM, remains unchanged.. To
compute CTM environments on the Kagome lattice, we coarse-grain three sites on
each down-pointing triangle into a single tensor resulting in an effective
square lattice. The local Hilbert space dimension thus grows to $3^{3}=27$,
which makes the optimizations memory-intensive. In Fig. 5, we demonstrate the
dramatic speed-up achieved by utilizing $U(1){\times}U(1)$ symmetry. For $D=9$
iPEPS, a single gradient step is already accelerated by more than a factor of
$100$. For larger bond dimensions, the simulations without symmetries become
prohibitive and we estimate the speed-up based on extrapolation of the scaling
at smaller bond dimensions.

Figure 5: Use case 2: optimization of $U(1){\times}U(1)$-symmetric iPEPS for
SU(3) Kagome model. In (a), scaling of the wall time per CTM step (in seconds)
for the entire gradient optimization step of iPEPS. The bond dimension of the
environment is $\chi=D^{2}$. The inset shows the relative speed-up compared to
an implementation without symmetry, with $D=12$ and $13$ simulation wall times
estimated from the extrapolation (blue dashed line). In (b), a distribution of
blocks of an enlarged corner by their size (number of elements) before the
fusion to a matrix. In (c), $L\times L$ blocks of the block-diagonal enlarged
corner after fusion. We plot them as a heatmap, with different $U(1)$ charges
on x- and y-axes.

In contrast to the example in Sec. III.1 utilizing the $U(1)$-symmetry, the
speed-up in this case is not monotonic in $D$. This happens due to the varying
structure of the iPEPS tensors, i.e., the allowed symmetry sectors and their
sizes. In Fig. 5(b,c) we illustrate the block structure of enlarged corners
before and after the fusion to a block-diagonal matrix. Generally, for larger
groups the number of blocks of enlarged corner before fusion is substantially
higher. Even at $D=7$, the total number of blocks is already more than $3000$
whereas for $U(1)$-symmetric enlarged corner in Sec. III.1 it was below $300$.
For $D=13$ ansatz, the fusion of enlarged corner into a block-diagonal matrix
requires processing of more than $16000$ blocks, with more than half of them
being small in size, having roughly $\mathcal{O}(10^{3})$ elements or less.
This granularity defines the bottleneck of the simulations. For $D=12$ and
$D=13$ the ratios between the computational time of SVD and contractions
including fusion to block-diagonal enlarged corners are 3:4 and 1:10
respectively. Overall, the $U(1){\times}U(1)$ simulations become dominated by
fusion, with SVD being subleading. The precise speed-up rests on the sizes of
blocks, such as here, where $D=12$ has a slightly higher proportion of largest
blocks compared to $D=13$.

Figure 6: Use case 3: expectation values from CTM environments in fermionic
iPEPS for finite-temperature Hubbard model. The bond dimension of the
environment is $\chi=5D$, which is sufficient here to converge the expectation
values. We show, in (a), the wall time per CTM step and, in (b), the data size
(memory requirement) of the biggest intermediate tensor appearing while
building enlarged CTM corners in Fig. 3(b). Panel (c) shows the number of
blocks in an enlarged corner before the fusion in Fig. 3(c), and panel (d) is
the size of the largest $L\times L$ block for SVD in Fig. 3(d).

###  III.3 2D Fermi-Hubbard model on a square lattice

We consider a two-dimensional Fermi-Hubbard model (FHM) with on-site repulsion
as studied in Ref. [48]. The Hamiltonian reads

$H=-t\sum_{\langle
i,j\rangle,\sigma}\left(c^{\dagger}_{i\sigma}c_{j\sigma}+\text{h.c.}\right)+U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\mu\sum_{i}n_{i},$
(17)

where $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ are the creation and
annihilation operators for an electron with spin $\sigma$ at site $i$,
$n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is a corresponding number
operator, $t$ is the hopping amplitude, $U$ is the on-site Coulomb repulsion
strength, and $\mu$ is the chemical potential.

The iPEPS ansatz employs a checkerboard lattice with a 2-site unit cell. The
thermal state for inverse temperature $\beta$ is obtained by evolving the
infinite-temperature purification $|\psi(0)\rangle$ under the propagator
$U(\beta)=e^{-\frac{\beta}{2}H}$. The initial purification
$|\psi(0)\rangle=\prod_{j}\prod_{m=\uparrow,\downarrow}\frac{1}{\sqrt{2}}\sum_{s_{m_{j}}=a_{m_{j}}=0,1}|s_{m_{j}}a_{m_{j}}\rangle$
is a product of maximally entangled pairs between each physical site and its
corresponding ancilla, translating to local Hilbert spaces of dimension
$4^{2}=16$. In the examples below, we run the imaginary time evolution
employing the NTU scheme targeting $\beta=2$.

In YASTN, the fermionic exchange order is implemented following the scheme of
Refs. [56, 57] by projecting the lattice ansatz onto a plane, imposing a
canonical fermionic order, and applying swap gates on every line crossing that
arises; see Fig. 3(b). Swap gate introduces sign changes for blocks with
charges of odd parity on both swapped legs. This makes $Z_{2}$ a minimal
symmetry needed for fermionic system simulations. The Hamiltonian in Eq. (17)
preserves the number of particles per spin direction, which allows us to
implement the model under $U(1){\times}U(1)$ symmetry as the highest abelian
symmetry.

The expectation values of the thermal state at final $\beta=2$ are calculated
using the CTM. Fig. 6 demonstrates an advantage of symmetric tensors by
comparing $Z_{2}$ (parity; minimal requirement for fermionic statistics),
$U(1)$ (total charge conservation), and $U(1){\times}U(1)$ (total charge
conservation for each spin) symmetries. We also show equivalent values for the
corresponding dense tensors with no symmetry. Here, we choose the
environmental bond dimension $\chi=5D$.

Fig. 6(a) presents the computational wall time for one CTM step as a function
of the iPEPS bond dimension for all the tested symmetries, with systematic
improvement offered by higher symmetries. Fig. 6(b) highlights the memory
usage bottleneck, showcasing the size of the largest object formed during the
CTM iteration, i.e., an intermediate step of the contraction in Fig. 3(b);
this tensor has to be later fused, and there are other tensors in the memory,
so the memory peak is roughly two times higher. It illustrates that those
simulations are ultimately memory-limited. Employing
$U(1){\times}U(1)$-symmetry offers a systematic $30$-fold memory gain as
compared to tensors with no symmetries, which ultimately allows for successful
simulations up to $D=36$.

Following the previous examples, in Fig. 6(c), we present the number of blocks
processed during the fusion that forms enlarged corners in Fig. 3(c). A
particular challenge for $U(1){\times}U(1)$ case is that the number of blocks
can exceed 10000. Nevertheless, the SVD is a dominant factor taking at least
half of the simulation time for $D\geq 25$ in our numerical experiments. In
Fig. 6(d), we show the (sectorial) bond dimension of the largest block
decomposed in Fig. 3(d), that is $\mathcal{O}(10^{4})$ for the largest bond
dimension achieved with each employed symmetry. The
$U(1){\times}U(1)$-symmetry offers here $15$-fold improvement as compared to a
setup with no symmetries involved.

##  IV Future outlook

Tensor networks are becoming increasingly popular tool for numerical treatment
of quantum systems, ranging from ground state simulations of condensed-matter
systems to simulation of quantum circuits. The landscape of associated
software is continuously growing. For 1D and quasi-1D geometries, well-
established and mature packages offer a rich set of MPS algorithms covering
direct energy minimization, (imaginary) time evolution, and much more. For
two-dimensional geometries, predominantly targeted by iPEPS, the field remains
nascent.

Here, we have introduced YASTN, a python-based TN library with strong emphasis
on simulation of two-dimensional systems by iPEPS, motivated by the need for
both abelian symmetries and automatic differentiation. By design, the dense
linear algebra and the AD engine are provided by different backends, allowing
for implementation-specific optimizations. YASTN, with its rich set of
examples covering ground state simulations of various 2D spin lattice models
(through peps-torch) and finite-temperature simulations of 2D Hubbard model,
thus joins similar efforts by VariPEPS [58], PEPSkit [59], ad-peps [36], and
peps-torch [35] together lowering the barrier for entry.

The wide separation between the high-level description of iPEPS algorithms and
their fast execution, optimized down to low-level dense linear algebra,
especially for symmetric tensors, remains a challenge. Unlike MPS simulations,
iPEPS contraction algorithms for computing environments and evaluation of
observables involve a more diverse set of tensor contractions, varying in
ranks and block sparsity patterns. Furthermore, flexible deployment and the
ability to leverage heterogenous clusters, accounting for the iPEPS-specific
block sparsity, is vital for addressing the sharp $\mathcal{O}(D^{8}\sim
D^{12})$ (albeit polynomial) scaling with the bond dimension, which is the key
resource governing the precision of iPEPS. These challenges thus call for
further development.

###### Acknowledgements.

We thank Philippe Corboz, Piotr Czarnik, Jacek Dziarmaga, Boris Ponsioen,
Yintai Zhang for inspiring discussions that were invaluable in the development
of this package. We acknowledge the funding by the National Science Center
(NCN), Poland, under projects 2019/35/B/ST3/01028 (A.S.), 2020/38/E/ST3/00150
(G.W.), and project 2021/03/Y/ST2/00184 within the QuantERA II Programme that
has received funding from the European Union Horizon 2020 research and
innovation programme under Grant Agreement No. 101017733 (M.M.R.). J.H.
acknowledges support from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 101001604) and from the Swiss National Science Foundation through
a Consolidator Grant (iTQC, TMCG-2_213805).

## References

  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. 349, 117 (2014).
  * Orús [2019] R. Orús, Tensor networks for complex quantum systems, Nat. Rev. Phys. 1, 538 (2019).
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * White [1993] S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993).
  * Rommer and Östlund [1997] S. Rommer and S. Östlund, Class of ansatz wave functions for one-dimensional spin systems and their relation to the density matrix renormalization group, Phys. Rev. B 55, 2164 (1997).
  * Dukelsky  _et al._ [1998] J. Dukelsky, M. A. Martín-Delgado, T. Nishino, and G. Sierra, Equivalence of the variational matrix product method and the density matrix renormalization group applied to spin chains, EPL 43, 457 (1998).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev. Lett. 93, 040502 (2004).
  * Schuch  _et al._ [2008] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Entropy scaling and simulability by matrix product states, Phys. Rev. Lett. 100, 030504 (2008).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions, arXiv:cond-mat/0407066 (2004).
  * Verstraete  _et al._ [2006] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006).
  * Orús [2012] R. Orús, Exploring corner transfer matrices and corner tensors for the classical simulation of quantum lattice systems, Phys. Rev. B 85, 205117 (2012).
  * Vlaar and Corboz [2021] P. C. G. Vlaar and P. Corboz, Simulation of three-dimensional quantum systems with projected entangled-pair states, Phys. Rev. B 103, 205137 (2021).
  * Vidal [2007] G. Vidal, Classical simulation of infinite-size quantum lattice systems in one spatial dimension, Phys. Rev. Lett. 98, 070201 (2007).
  * Jordan  _et al._ [2008] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Phys. Rev. Lett. 101, 250602 (2008).
  * McCulloch [2007] I. P. McCulloch, From density-matrix renormalization group to matrix product states, J. Stat. Mech. Theory Exp. 2007, P10014 (2007).
  * Singh  _et al._ [2010] S. Singh, R. N. C. Pfeifer, and G. Vidal, Tensor network decompositions in the presence of a global symmetry, Phys. Rev. A 82, 050301 (2010).
  * Singh  _et al._ [2011] S. Singh, R. N. C. Pfeifer, and G. Vidal, Tensor network states and algorithms in the presence of a global U(1) symmetry, Phys. Rev. B 83, 115125 (2011).
  * Singh and Vidal [2012] S. Singh and G. Vidal, Tensor network states and algorithms in the presence of a global SU(2) symmetry, Phys. Rev. B 86, 195114 (2012).
  * Silvi  _et al._ [2019] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes , 8 (2019).
  * [22] M. M. Rams, G. Wójtowicz, A. Sinha, and J. Hasik, Yastn: Yet another symmetric tensor network, https://github.com/yastn/yastn.
  * Corboz [2016] P. Corboz, Variational optimization with infinite projected entangled-pair states, Phys. Rev. B 94, 035133 (2016).
  * Vanderstraeten  _et al._ [2016] L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Gradient methods for variational optimization of projected entangled-pair states, Phys. Rev. B 94, 155123 (2016).
  * Liao  _et al._ [2019] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Differentiable programming tensor networks, Phys. Rev. X 9, 031041 (2019).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), SciPost Phys. Lect. Notes , 5 (2018).
  * Zhai  _et al._ [2023] H. Zhai, H. R. Larsson, S. Lee, Z.-H. Cui, T. Zhu, C. Sun, L. Peng, R. Peng, K. Liao, J. Tölle, J. Yang, S. Li, and G. K.-L. Chan, Block2: A comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond, J. Chem. Phys. 159, 234801 (2023).
  * Ballarin  _et al._ [2024] M. Ballarin, G. Cataldi, A. Costantini, D. Jaschke, G. Magnifico, S. Montangero, S. Notarnicola, A. Pagano, L. Pavesic, M. Rigobello, N. Reinić, S. Scarlatella, and P. Silvi, Quantum tea: qtealeaves (2024).
  * [30] C. R. Roberts _et al._ , Tensornetwork: A library for easy and efficient manipulation of tensor networks, https://github.com/google/TensorNetwork.
  * Wu  _et al._ [2023] K.-H. Wu, C.-T. Lin, K. Hsu, H.-T. Hung, M. Schneider, C.-M. Chung, Y.-J. Kao, and P. Chen, The cytnx library for tensor networks, SciPost Phys. Codebases (2023).
  * Motoyama  _et al._ [2022] Y. Motoyama, T. Okubo, K. Yoshimi, S. Morita, T. Kato, and N. Kawashima, Tenes: Tensor network solver for quantum lattice systems, Comput. Phys. Commun. 279, 108437 (2022).
  * Haegeman and contributors [2024] J. Haegeman and contributors, TensorKit.jl (2024).
  * Weichselbaum [2024] A. Weichselbaum, Qspace - an open-source tensor library for Abelian and non-Abelian symmetries, arXiv:2405.06632 (2024).
  * [35] J. Hasik and contributors, peps-torch: Solving two-dimensional spin models with tensor networks (powered by pytorch), https://github.com/jurajHasik/peps-torch.
  * [36] B. Ponsioen, ad-peps: ipeps ground- and excited-state implementation based on automatic differentiation, https://github.com/b1592/ad-peps.
  * Naumann  _et al._ [2024] J. Naumann, E. L. Weerda, M. Rizzi, J. Eisert, and P. Schmoll, An introduction to infinite projected entangled-pair state methods for variational ground state simulations using automatic differentiation, arXiv:2308.12358 (2024).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, _et al._ , Array programming with numpy, Nature 585, 357 (2020).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , Pytorch: An imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Pfeifer  _et al._ [2014] R. N. C. Pfeifer, G. Evenbly, S. Singh, and G. Vidal, Ncon: A tensor network contractor for matlab, arXiv:1402.0939 (2014).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * F. Verstraete and Cirac [2008] V. M. F. Verstraete and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).
  * King  _et al._ [2024] A. D. King, A. Nocera, M. M. Rams, J. Dziarmaga, R. Wiersema, W. Bernoudy, J. Raymond, N. Kaushal, N. Heinsdorf, R. Harris, _et al._ , Computational supremacy in quantum simulation, arXiv:2403.00910 (2024).
  * Wójtowicz  _et al._ [2023] G. Wójtowicz, A. Purkayastha, M. Zwolak, and M. M. Rams, Accumulative reservoir construction: Bridging continuously relaxed and periodically refreshed extended reservoirs, Phys. Rev. B 107, 035150 (2023).
  * Dziarmaga [2021] J. Dziarmaga, Time evolution of an infinite projected entangled pair state: Neighborhood tensor update, Phys. Rev. B 104, 094411 (2021).
  * Dziarmaga [2022] J. Dziarmaga, Simulation of many-body localization and time crystals in two dimensions with the neighborhood tensor update, Phys. Rev. B 105, 054203 (2022).
  * Czarnik  _et al._ [2019] P. Czarnik, J. Dziarmaga, and P. Corboz, Time evolution of an infinite projected entangled pair state: An efficient algorithm, Phys. Rev. B 99, 035115 (2019).
  * Sinha  _et al._ [2022] A. Sinha, M. M. Rams, P. Czarnik, and J. Dziarmaga, Finite-temperature tensor network study of the hubbard model on an infinite square lattice, Phys. Rev. B 106, 195105 (2022).
  * Sinha  _et al._ [2024] A. Sinha, M. M. Rams, and J. Dziarmaga, Efficient representation of minimally entangled typical thermal states in two dimensions via projected entangled pair states, Phys. Rev. B 109, 045136 (2024).
  * Nishino and Okunishi [1996] T. Nishino and K. Okunishi, Corner transfer matrix renormalization group method, J. Phys. Soc. Japan 65, 891 (1996).
  * Orús and Vidal [2009] R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009).
  * Corboz  _et al._ [2014] P. Corboz, T. M. Rice, and M. Troyer, Competing states in the $t$-$j$ model: Uniform $d$-wave state versus stripe state, Phys. Rev. Lett. 113, 046402 (2014).
  * Fishman  _et al._ [2018] M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for contracting infinite two-dimensional tensor networks, Phys. Rev. B 98, 235148 (2018).
  * Hasik  _et al._ [2022] J. Hasik, G. B. Mbeng, S. Capponi, F. Becca, and A. M. Läuchli, Symmetric projected entangled-pair states analysis of a phase transition in coupled spin-$\frac{1}{2}$ ladders, Phys. Rev. B 106, 125154 (2022).
  * Xu  _et al._ [2023] Y. Xu, S. Capponi, J.-Y. Chen, L. Vanderstraeten, J. Hasik, A. H. Nevidomskyy, M. Mambrini, K. Penc, and D. Poilblanc, Phase diagram of the chiral SU(3) antiferromagnet on the kagome lattice, Phys. Rev. B 108, 195153 (2023).
  * Corboz and Vidal [2009] P. Corboz and G. Vidal, Fermionic multiscale entanglement renormalization ansatz, Phys. Rev. B 80, 165129 (2009).
  * Corboz  _et al._ [2010] P. Corboz, R. Orús, B. Bauer, and G. Vidal, Simulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair states, Phys. Rev. B 81, 165104 (2010).
  * Naumann  _et al._ [2023] J. Naumann, E. L. Weerda, M. Rizzi, J. Eisert, and P. Schmoll, varipeps – a versatile tensor network library for variational ground state simulations in two spatial dimensions, arXiv:2308.12358 (2023).
  * [59] M. van Damme and contributors, PEPSKit.jl: Tools for working with projected entangled-pair states, https://github.com/quantumghent/PEPSKit.jl.


# Fabrication of thin diamond membranes by Ne+ implantation

Luca Basso<EMAIL_ADDRESS>Center For Integrated Nanotechnology, Sandia
National Laboratories, Albuquerque, New Mexico 87123, USA  Michael Titze
Sandia National Laboratories, Albuquerque, New Mexico 87185, USA  Jacob
Henshaw  Center For Integrated Nanotechnology, Sandia National Laboratories,
Albuquerque, New Mexico 87123, USA  Pauli Kehayias  Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  Rong Cong  Department of
Physics, Brown University, Providence, Rhode Island 02912, USA  Maziar Saleh
Ziabari  Center for High Technology Materials and Department of Physics and
Astronomy, University of New Mexico, Albuquerque, New Mexico 87131, USA  Tzu-
Ming Lu  Center For Integrated Nanotechnology, Sandia National Laboratories,
Albuquerque, New Mexico 87123, USA  Michael P. Lilly  Center For Integrated
Nanotechnology, Sandia National Laboratories, Albuquerque, New Mexico 87123,
USA  Andrew M. Mounce<EMAIL_ADDRESS>Center For Integrated
Nanotechnology, Sandia National Laboratories, Albuquerque, New Mexico 87123,
USA

###### Abstract

Color centers in diamond are one of the most promising tools for quantum
information science. Of particular interest is the use of single-crystal
diamond membranes with nanoscale-thickness as hosts for color centers. Indeed,
such structures guarantee a better integration with a variety of other quantum
materials or devices, which can aid the development of diamond-based quantum
technologies, from nanophotonics to quantum sensing. A common approach for
membrane production is what is known as “smart-cut”, a process where membranes
are exfoliated from a diamond substrate after the creation of a thin sub-
surface amorphous carbon layer by He+ implantation. Due to the high ion
fluence required, this process can be time-consuming. In this work, we
demonstrated the production of thin diamond membranes by neon implantation of
diamond substrates. With the target of obtaining membranes of $\sim$200 nm
thickness and finding the critical damage threshold, we implanted different
diamonds with 300 keV Ne+ ions at different fluences. We characterized the
structural properties of the implanted diamonds and the resulting membranes
through SEM, Raman spectroscopy, and photoluminescence spectroscopy. We also
found that a SRIM model based on a two-layer diamond/sp2-carbon target better
describes ion implantation, allowing us to estimate the diamond critical
damage threshold for Ne+ implantation. Compared to He+ smart-cut, the use of a
heavier ion like Ne+ results in a ten-fold decrease in the ion fluence
required to obtain diamond membranes and allows to obtain shallower smart-
cuts, i.e. thinner membranes, at the same ion energy.

##  I Introduction

Diamond is one of the leading material platforms for quantum information
science (QIS) not only for its unique material properties, but for the
presence of optically-active lattice defects, known as color centers. Indeed,
these diamond defects’ properties - such as spin-dependent optical transitions
[1, 2], exceptional coherence times [3, 4], single photon emission [5, 6], and
controllable interaction with the surrounding environment [7, 8, 9] \- have
made them ideal candidates for different QIS applications like quantum
computing [10, 11], networking [12], and sensing [13, 14, 15, 16]. Most of
these applications can benefit by the use of color centers embedded in a sub-
micrometer thick, high-quality diamond membrane. For instance, these membranes
can lead to an improved coupling of single emitters to waveguides or cavities
[17, 18, 19], an increased entanglement efficiency [20], and a better
interface with other materials for quantum sensing [21, 22]. In general,
color-center-enriched diamond membranes can advance the development of
diamond-based QIS technologies thanks to an improved integration with a
variety of quantum materials or devices [23]. One of the most significant
hurdles for the development of diamond-membrane-based quantum technologies is
membrane production, since growth of high crystal quality, thin diamond film
on non-diamond substrate is still challenging [24, 25]. For this reason, the
most common technique for membrane production is diamond “smart-cut” used in
combination with high-quality chemical vapor deposition (CVD) diamond
overgrowth [23, 26, 27, 28, 29]. Smart-cut involves implanting a bulk diamond
plate with He+ ions to create an amorphous layer beneath diamond surface.
Following this, a single-crystal, thin diamond film is grown over the
implanted diamond. The amorphous layer is then selectively etched resulting in
the exfoliation of the diamond membrane. Incorporation of color centers in the
overgrown layer, either by doping during growth or subsequent ion
implantation, has already been demonstrated for NV centers [30] as well as
other group IV defects such as SiV [31], GeV [23] and SnV [32] centers with
bulk-like spin coherence properties. However, due to the low mass of He,
diamond smart-cut requires extremely high ion fluence, in the range
$(0.5-1.5)\times 10^{17}$ ions/cm2 [23, 33].

To make this process more efficient, by reducing beam-time and cost, the
smart-cut can be performed with a heavier ion, such as neon. In this work, we
report the fabrication of ultrathin diamond membranes through Ne+ implantation
of diamond. We used Ne+ instead of other heavier noble gases because smart-cut
requires substrate atoms and irradiating ions to have comparable size. Larger
ions, such as Ar+, Kr+ or Xe+, do not penetrate the surface easily, leading to
structural damage limited at the target surface resulting in an increased
sputtering yield [34]. Moreover, Ne+ does not introduce unwanted color centers
in diamond, as Xe does [35]. To obtain $\sim$200 nm membranes we implanted
diamonds with different ion beam fluences in the range $(1.38-11.0)\times
10^{15}$ ions/cm2 with 300 keV energy. The implanted diamonds were
characterized through SEM to observe the depth of the amorphous layer, and
with Raman and photoluminescence (PL) spectroscopy to study the structural
properties of the thin diamond layer capping the amorphous region. The diamond
membranes were then isolated by removing the amorphous carbon layer through
electrochemical (EC) etching in a boric acid solution. Finally, we estimated
the damage threshold $D_{c}$, i.e. the minimum ion-induced vacancy density
required to induce diamond amorphization, by addressing the discrepancies
between the experimental results with Stopping Range of Ions in Matter (SRIM)
simulations. We found that a target modeled as a two layers diamond/sp2-carbon
system better describes ion implantation with SRIM. From the obtained $D_{c}$
value, we estimated the lowest fluence required to obtain a smart-cut is
$(2.3\pm 0.1)\times 10^{15}$ ions/cm2 (the lowest fluence smart-cut
experimentally demonstrated in this work is of $2.75\times 10^{15}$ ions/cm2),
more than an order of magnitude smaller compared to He+ smart-cut. A faster
and more efficient diamond membrane smart-cut production process could benefit
and drive the development of a variety of diamond-based quantum technologies.

Figure 1: (a) Schematics of Ne+ smart-cut process. (1) Ne+ implantation to form a damaged layer (dark grey layer) buried below diamond surface. (2) High temperature annealing to convert the damaged layer into an etchable amorphous sp2 carbon phase (black layer). (3) Electrochemical etching in boric acid solution to etch the sp2 carbon layer. (4) Membranes are partially exposed, i.e. the amorphous layer is removed only in a section of the sample, and still attached to the bulk diamond. (b) SRIM simulation of vacancy density profiles for a 300 keV Ne+ implantation with different fluences. The critical damage threshold range $D_{c}=1-9\times 10^{22}$ vacancies/cm3 for diamond He+ implantation estimated from the literature is indicated by the grey box (${D_{c}}$ ref). The green line (${D_{c}}$ this work) is the value estimated in this paper of $D_{c}=(8.4\pm 0.8)\times 10^{22}$ vacancies/cm3. Sample |  Ne+ fluence |  Amorphous |  Top diamond |  Diamond
---|---|---|---|---
|  ($10^{15}$ ions/cm2) |  layer (nm) |  layer before EC etching (nm) |  membrane (nm)
1 |  1.38 |  |  |
2 |  2.75 |  $262\pm 12$ |  $219\pm 12$ |  $222\pm 11$
3 |  4.12 |  $338\pm 12$ |  $203\pm 12$ |  $195\pm 10$
4 |  5.50 |  $363\pm 12$ |  $169\pm 12$ |  $176\pm 11$
5 |  6.88 |  $613\pm 15$ |  |  $61\pm 8$
6 |  8.25 |  $670\pm 15$ |  |  $55\pm 6$
7 |  9.62 |  $676\pm 15$ |  |  $39\pm 4$
8 |  11.0 |  $741\pm 16$ |  |
Table 1: List of samples with the relative Ne+ implantation fluence and
details of SEM analysis before and after EC etching. For sample 1, the
amorphous layer is not formed due to a low ion fluence. For samples 2-8 the
diamond layer above the amorphous region, i.e. “Diamond layer before EC
etching” could not be observed with SEM imaging. For sample 8 the
amorphization reached the surface and no membrane is observed after the EC
etching.

##  II Experimental Methods

###  II.1 Sample preparation

####  II.1.1 Ion implantation and thermal annealing

Diamond membranes were obtained by implantation of 3$\times$3$\times$0.25 mm3
optical-grade ($\leq 1$ ppm [N]) single-crystal diamonds produced by Element
Six. Diamond plates were implanted with 20Ne+ with an energy of 300 keV with a
fluence ranging from 1.38$\times 10^{15}$ to 11$\times 10^{15}$ ions/cm2 (see
Table 1 for the list of all fluences). Implantation is performed using a 3 MV
Pelletron accelerator using a plasma source containing a mixture of H and Ne.
Samples are held perpendicular to the incident ion beam by means of carbon
tape. The sample temperature is maintained $<25^{\circ}$C through cooling by a
liquid nitrogen cold-finger cryostat to avoid melting of the carbon tape from
the high power of the ion beam. The ion beam is scanned over a 1/4” (6.35 mm)
mask to ensure a homogeneous implantation profile. The fluence is measured
through current integration off the copper sample holder. The samples were
then annealed at 1200 ∘C for 2h (heating rate of 2.5 ∘C/min) in ultra-high
vacuum.

####  II.1.2 Electrochemical (EC) etching

The implanted diamonds were glued in the center of a Petri dish, which was
then filled with a saturated solution of boric acid in deionized (DI) water
[36]. Two copper electrodes with surface area $\sim$ 1 cm2 were immersed in
the solution and brought as close as possible to the two opposite sides of the
diamond, to a distance between the anode and the cathode of $\lesssim$ 1 cm.
We then used a DC power supply to apply 300 V between the electrodes. The
concentration of the solution was adjusted by adding DI water or boric acid
solution to maintain a current across the sample in the range $15-20$ mA. The
samples were only partially etched, resulting in the membranes still being
attached to the bulk diamond after the partial removal of the amorphous layer.
After the etching, the diamonds were rinsed in acetone and then cleaned in a
piranha solution to remove any contaminants from the EC etching.

Figure 2: (a) Schematics of the SEM imaging geometry. The samples are titled
by a 52∘ angle, so that the surface is perpendicular respect to the FIB (green
beam) that is used to mill the implanted diamonds and expose the amorphous
layer (black layer). SEM (red beam) is then used to image the tilted samples
so the cross section can be observed. (b) SEM images of the smart-cuts
obtained with different Ne+ fluences (reported on top of each picture). White
dashed lines indicate the amorphous layer position. Scale bar of 500 nm is the
same for all the pictures. (c) Amorphous carbon layer thickness (green points)
and thickness of the top diamond layer (blue points) measured from (b) as
function of Ne+ fluence. Values used in this plot are reported in Table 1.

####  II.1.3 SRIM simulation

To simulate the defect density profile generated by Ne+ implantation of
diamond, we used Stopping Range of Ions in Matter (SRIM) Monte Carlo
simulation software [37]. To obtain the vacancy density profile we simulated
10,000 Ne+ ions with 19.992 amu mass, 300 keV energy, and 500 nm target depth.
The resulting vacancy density curve was then multiplied by the different ion
fluences to obtain the different vacancy density profiles as function of Ne+
implantation fluence. The diamond target was modeled as a 12C solid with: 3.51
g/cm3 density, 37.5 eV atom displacement threshold energy, 3 eV lattice damage
threshold, and 7.4 eV surface damage threshold [38]. To model the two-layer
diamond/sp2-carbon target with SRIM, we set the diamond with the same
parameters as above, while the sp2-carbon ones were changed into the standard
SRIM parameters for graphite: 2.253 g/cm3 density, 28 eV atom displacement
threshold energy, 3 eV lattice damage threshold, and 7.41 eV surface damage
threshold.

###  II.2 Sample characterization

####  II.2.1 SEM and FIB

FIB milling and SEM imaging were performed with a FEI Nanolab 650 microscope.
We used a Ga+ FIB accelerated at 30 kV with a current of 3 nA to obtain a
substantial diamond etching rate, whereas SEM was operated in secondary
electron mode with an acceleration voltage of 5 kV. To reduce charging effects
during the analysis, a layer of Pt was deposited over the region of interest
before FIB milling. This layer also helped to clearly determine the position
of the sample surface during SEM imaging of the diamond cross section, thus
allowing to correctly determine the damaged layer depth.

####  II.2.2 Raman spectroscopy

Raman and PL analysis were obtained with a WITec alpha300 Raman system. The
excitation was provided by a 532 nm laser with 3.4 mW output power which was
focused on the sample surface with a 100$\times$ 0.9 NA objective to a laser
spot diameter of $\sim$ 20 $\upmu$m. The emitted light was dispersed with a
1800 g/mm grating and collected by a CCD camera.

##  III Results and discussion

###  III.1 Ne smart-cut

Schematics of the Ne+ smart-cut process are given in Fig. 1(a). The diamond
substrates were implanted with Ne+ ions to form a damaged layer buried below
diamond surface. We then annealed the samples to convert the region damaged
above the critical damage threshold $D_{c}$ into an etchable sp2 amorphous
carbon layer. The substrates were finally electrochemically etched to
selectively remove the amorphous carbon layer and release the diamond
membrane. To obtain different membrane thicknesses, we smart-cut nine diamonds
with different Ne+ fluences at the same energy of 300 keV. We set the
different fluences, ranging from 1.38$\times 10^{15}$ to 11$\times 10^{15}$
ions/cm2, by comparing the SRIM simulation [37] for Ne+ implantation of
diamond with the critical damage threshold range for He+ implantation reported
in the literature $D_{c}=(1-9)\times 10^{22}$ vacancies/cm3 [39], as shown in
Fig. 1(b). The samples, labeled from 1 to 9, and the relative Ne+ implantation
fluence are reported in Table 1.

Figure 3: (a) Raman spectra of implanted diamonds with different Ne+ fluences.
Solid lines are the Lorentzian fit of the experimental data (circles) with a
linear background. Each curve is normalized respect to the maximum peak value.
The black vertical dashed line is the position of the Raman peak for
unimplanted diamond. (b) Comparison between the diamond Raman peak of sample 1
(dark green) and unimplanted diamond (dark red), both fit with a Lorentzian
function. The Raman line for unimplanted diamond is centered at $1330.1\pm
0.1$ cm-1 with a FWHM of $2.7\pm 0.1$ cm-1, whereas sample 1 shows the Raman
peak at $1329.9\pm 0.1$ cm-1 with a FWHM of $2.8\pm 0.1$ cm-1. (c) Raman
spectrum of sample 9. Data points (light blue circles) are fit by a five peaks
model [40] with the peaks labelled as: G at $\sim$ 1600 cm-1, D1 $\sim$ 1346
cm-1, D2 $\sim$ 1648 cm-1, D3 $\sim$ 1532 cm-1, D4 $\sim$ 1202 cm-1. (d) PL
spectra from native NV centers in the Ne+ implanted diamonds compared to the
unimplanted one (grey curve).

###  III.2 Implanted diamonds analysis

To measure the thickness and the depth of the amorphous layer, we FIB-milled
the implanted substrate until the damaged layer was exposed, then used SEM to
image the diamond cross section, as shown in Fig. 2(a). The SEM images of the
implanted and annealed diamonds are reported in Fig. 2(b). Sample 1, implanted
with the lowest fluence of 1.38$\times 10^{15}$ ions/cm2, does not show the
presence of a damaged layer, meaning that the damage threshold $D_{c}$ was not
reached. In the other implanted samples the amorphous layers are visible as a
darker layer [39], with thicknesses increasing with Ne+ fluence from $262\pm
12$ nm to $741\pm 16$ nm, as shown in Fig. 2(c). For samples 2, 3 and 4, the
top diamond layer thickness above the amorphous one ranges from $219\pm 12$ nm
to $169\pm 12$ nm (see Table 1 for all measured values). For the remaining
samples 5 to 8, the top diamond layer cannot be observed. This is likely a
consequence of its very small thickness and the difficulty to SEM image small
features in insulating materials.

Figure 4: (a) SEM images of the resulting diamond membranes obtained with
different Ne+ fluence (reported on top of each picture). (b) Membrane
thickness (blue points) measured from (a) as function of Ne+ fluence.

To determine if there is a diamond layer above the amorphous carbon layer, we
used Raman spectroscopy and photoluminescence (PL) from the native NV centers.
The Raman spectra are shown in Fig. 3(a): for each spectrum the fit (solid
curves) of the experimental data is a Lorentzian peak with a linear
background, while the vertical dashed black line is the peak position of
unimplanted diamond at $1330.1\pm 0.1$ cm-1 (from Fig. 3(b)). While the sample
1 peak position is consistent with the peak position of unimplanted diamond,
the other samples show a shift of the peak position towards smaller
wavenumbers as a consequence of tensile stress acting on the diamond layer
[41]. This stress may originate from the swelling of the implanted layer
resulting from a density decrease after ion implantation [42]. Indeed, diamond
density reduction in proximity of the amorphous layer was previously observed
in He+ smart-cut diamond [39] and, even if this decreased density is not
enough to cause amorphization, it can cause lattice distortion and strain. The
width of the peak also increases with Ne+ fluence, as a result of increasing
damage in the crystal structure [43]. Moreover, with increasing fluence the
peak intensity decreases, as the thickness of the diamond layer where the
signal originates decreases. Samples 5, 6, and 7, for which a diamond layer
above the amorphous one could not be observed with SEM, still show a diamond
Raman peak (purple, pink, and light green data respectively), demonstrating
the presence of a diamond layer at the sample surface. Sample 8 (light blue
curve) does not show the presence of a peak that can be attributed to diamond
phase.

Figure 5: (a) Raman spectra of the resulting diamond membranes. Experimental
data (circles) are fit with two Lorentzian, one for the signal from the
diamond substrate (grey curve at $\sim$ 1330 cm-1) the other for the membrane
Raman signal (black curve at $\sim$ 1315 cm-1). The latter peaks parameters
are reported as a function of Ne+ fluence in (b) for peak intensity, (c) peak
position, and (d) peak width. The three black dashed lines represent the
values for unimplanted diamond with a Raman peak at $1330.1\pm 0.1$ cm-1
having FWHM of $2.7\pm 0.1$ cm-1. In picture (b) the peaks intensities for
each sample are normalized respect to the peak intensity of unimplanted
diamond. For (b-d) error bars are not visible when smaller compared to the
data point size.

To confirm that there is no damage-induced phase transition occurring in
sample 1, Fig. 3(b) shows its Raman spectrum (dark green curve) with an
unimplanted diamond (dark red curve). The two peaks are completely
overlapping, confirming the absence of damage and sp2 carbon phase in sample
1. Fig. 3(c) reports the Raman spectrum on a wider range (1000-1800 cm-1) for
sample 8. The fit is a result of a five-Gaussian model, typical of amorphous
carbon [40], where the five peaks (black curves labelled G, D1, D2, D3, and
D4) originate from different sp2 carbon atoms vibrational modes [44]. The lack
of a peak attributable to diamond in Fig. 3(c) demonstrates the absence of a
diamond layer on top of the damage-induced amorphous region.

Fig. 3(d) shows the PL spectra from the native NV centers of the implanted
diamonds compared with the NV emission from an unimplanted diamond surface
(grey curve). Sample 1 (dark green curve) still shows the zero phonon lines of
NV0 and NV- at respectively 575 nm and 637 nm, whereas the maximum of the
characteristic phonon sideband, usually at $\sim$ 680 nm, is red-shifted to
$\sim$ 750 nm. This can be attributed to the presence of neutral vacancy (GR1)
defects in the implanted diamond [45, 46]. The sharp peak at 572.5 nm is the
Raman signal of diamond (corresponding to 1330 cm-1 at 532 nm laser
excitation). The NV PL emission spectra from samples 2 to 8 are notably
different. The PL emission from samples 2-4 is consistent with the NV0 with
phonon sideband maximum at $\sim$ 650 nm. The suppression of the NV- emission
is consistent with the NVs being in proximity of sp2 carbon atoms and is often
observed from NVs in nanodiamonds [47, 48]. The reason is the introduction of
acceptor states by the sp2 defects causing Fermi-level pinning below the
NV-/NV0 transition [49], making the negative charge state of the NV unstable
and the NV population dominated by NV0. Furthermore, the PL intensity
decreases with Ne+ fluence as the active material, i.e. the diamond layer
above the amorphous region, gets thinner, until a PL is detected but the NV
features can no longer be identified, as observed for samples 5-9. Samples 2
to 9 also show the Raman features from the amorphous layer, the two broad
peaks in the region $550-600$ nm ($615-2130$ cm-1), superimposed over the
sharp Raman peak of diamond at 572.5 nm. As already observed in the Raman
spectra, the intensity of the diamond Raman peaks decreases with Ne+ fluence.
The absence of the Raman signal from amorphous carbon in the PL spectrum of
sample 1 is an additional proof that the 1.38$\times 10^{15}$ ions/cm2
implantation does not cause a phase transition and the increased vacancy
concentration observed in the PL phonon sideband is not enough to damage the
lattice.

###  III.3 Diamond membrane characterization

To release the diamond membrane, we removed the amorphous carbon layer by EC
etching. The samples were partially etched, so that the membranes were still
attached on the diamond substrates after the EC etching, as shown in Fig.
1(a). SEM images of the resulting membranes cross sections are shown in Fig.
4(a). The membranes thicknesses decrease with Ne+ implantation fluence, from
$222\pm 11$ nm for sample 2, to $39$ $\pm$ $4$ nm for sample 8. Membrane
thickness as function of Ne+ fluence is reported in Fig. 4(b). Moreover, we
could also observe that the membrane thickness measured before the EC etching,
blue data in Fig. 2(c), matches the thickness of the final membranes for
sample 2, 3 and 4.

Results for the structural analysis performed with Raman spectroscopy are
shown in Fig. 5(a). The Raman spectra are dominated by the signal coming from
the diamond substrate, fit with a Lorentzian (grey curve) at $\sim$ 1330 cm-1,
as after the EC etching the membranes are laying on the bulk diamond
substrate. The Raman signal originating from the membranes is the smaller peak
at lower wavenumbers that is fit with a Lorentzian (black curve at $\sim$ 1315
cm-1). Details of intensity, position, and width of these Raman peaks are
plotted as function of Ne+ fluence in Fig. 5(b), (c), and (d) respectively,
together with the reference value of unimplanted diamond (dashed black lines)
from Fig. 3(b). From Fig. 5(b) we could observe that the membrane Raman peak
decreases with Ne+ fluence, as expected due to a decreasing material
thickness. Fig. 5(c) shows that the tensile stress is still present on the
membranes, due to a shift of peak position towards smaller wavenumbers, as
already observed in Fig. 3(a). We could observe that the shift increases with
Ne+ fluence, as for thinner membranes the strained section represent a larger
portion of the sample, thus leading to an increasing shift in the Raman peak.
Furthermore, the width of the membranes Raman peaks is $\sim 10\times$ larger
compared to the FWHM of the unimplanted diamond Raman peak, a signal of
partially damaged crystal structures. Shifts and broadening of similar
amplitudes have also been observed for Raman signal of He+ smart-cut diamonds
membranes [28, 50].

Figure 6: (a) SEM images of diamond implanted with 300 keV angle at
$5.50\times 10^{15}$ ions/cm2 with a 7∘ angle (top) and of the resulting
membrane after EC etching (bottom). White dashed lines indicate the amorphous
layer position. (b) Comparison between the measured amorphous layer thickness
and depth (grey area) with two SRIM simulations: one for Ne+ in diamond (red
curve), one for Ne+ implantation of a two layers target made of 113 nm of
diamond and 500 nm graphite (blue curve). Green area is the estimated damage
critical threshold of $D_{c}=(8.4\pm 0.8)\times 10^{22}$ vacancies/cm3.

###  III.4 Critical damage threshold estimation

When observing the amorphous carbon layer thickenss and depth, we could
observe that simulations (Fig. 1(b)) and experimental results (Fig. 2(a)) do
not match, in particular with SRIM underestimating the depth of the damaged
layers. As a matter of fact, SRIM does not take into account some phenomena
that may occur during ion implantation, such as channeling [51] or the change
of material properties. To check whether channeling plays a relevant role, we
implanted another diamond tilted by 7∘, with the same energy of 300 keV and a
Ne+ fluence of 5.5$\times 10^{15}$ ions/cm2. The SEM cross-sectional images
are reported in Fig. 6(a), after the implantation (top) and after the
amorphous layer is etched (bottom). The amorphous layer is $282\pm 13$ nm,
capped by a diamond layer of $113\pm 10$ nm; after the EC etching the
thickness of the resulting membrane is $121\pm 9$ nm. Fig. 6(b) compares this
experimental result (grey area) with SRIM simulation of Ne+ in diamond (red
curve). Despite the 7∘ implantation angle, the simulation still underestimates
the amorphous layer thickness, meaning that channeling does not play a major
role. To further asses this discrepancy we performed a SRIM simulation, blue
curve in Fig. 6(b), of Ne+ implanted in a two-layer substrate made of 113 nm
of diamond and 500 nm of sp2-carbon with same energy and fluence. This two-
layer target is a similar composition of the one expected at the end of the
implantation process, and since sp2-carbon has a lower density than diamond,
this simulation should overestimate the damaged layer depth. This was actually
observed, as shown in Fig. 6(b), which shows that the trailing edge of the
amorphous layer lays in between the two SRIM simulation tails. Thus, SRIM
inaccuracy in predicting the amorphous layer depth and thickness may be a
consequence of the change in material physical properties and phase transition
occurring during the implantation. Lastly, from comparing the SEM images and
the SRIM simulations we could estimate the critical damage threshold $D_{c}$,
resulting in $D_{c}=(8.4\pm 0.8)\times 10^{22}$ vacancies/cm3 (green area in
Fig. 6(b)). This value is consistent with the literature reported results for
$D_{c}$ in diamond, that for He+ implantation is in the range $1-9\times
10^{22}$ vacancies/cm3 [39]. From this value of $D_{c}$ we used SRIM to
estimate the minimum fluence $F_{min}$ for Ne+ at 300 keV required to obtain a
smart-cut (defined as the Ne+ fluence resulting in a damage profile where the
maximum value is equal to $D_{c}$), obtaining $F_{min}=(2.3\pm 0.1)\times
10^{15}$ ions/cm2.

##  IV Conclusions

In summary, we have demonstrated the fabrication of thin diamond membranes by
Ne+ implantation of diamond plates. The resulting membranes have tunable
thicknesses, ranging from $39\pm 4$ nm to $222\pm 11$ nm, depending on the Ne+
ions fluence. Moreover, by addressing the discrepancy between the experimental
results and SRIM simulation, we obtained an estimate for diamond amorphization
threshold. The use of a heavier ion, compared to the standard He+ “smart-cut”,
has two main advantages. Firstly, it allows to reduce the implantation fluence
required to form the amorphous layer by more then a factor 10. Furthermore, at
the same ion energy Ne+ has a smaller implantation depth compared to He,
leading to a shallower damaged layer in turn resulting in a thinner diamond
membrane. These two points, when considering that these smart-cut membranes
are a sacrificial layer etched after being used as seeds for high-quality
diamond overgrowth, show the relevance of our findings. A faster and more
efficient method to produce thin diamond membranes to use as CVD templates for
high-quality diamond overgrowth may benefit a wide range of diamond membranes
applications, such as quantum information processing, nanophotonics, or
diamond-based electronics.

## Acknowledgements

We thank Ania C. Bleszynski Jayich, Dirk Bouwmeester, Hyunseok Oh, and Ian
Hedgepeth from UCSB for fruitful discussion. Sandia National Laboratories is a
multi-mission laboratory managed and operated by National Technology and
Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell
International, Inc., for the DOE’s National Nuclear Security Administration
under contract DE-NA0003525. This work was funded, in part, by the Laboratory
Directed Research and Development Program and performed, in part, at the
Center for Integrated Nanotechnologies, an Office of Science User Facility
operated for the U.S. Department of Energy (DOE) Office of Science. This paper
describes objective technical results and analysis. Any subjective views or
opinions that might be expressed in the paper do not necessarily represent the
views of the U.S. Department of Energy or the United States Government.

## References

## References

  * [1] Doherty M W, Manson N B, Delaney P, Jelezko F, Wrachtrup J and Hollenberg L C 2013 Phy. Rep. 528 1–45
  * [2] Rogers L J, Jahnke K D, Metsch M H, Sipahigil A, Binder J M, Teraji T, Sumiya H, Isoya J, Lukin M D, Hemmer P et al. 2014 Phys. Rev. Lett. 113 263602
  * [3] Balasubramanian G, Neumann P, Twitchen D, Markham M, Kolesov R, Mizuochi N, Isoya J, Achard J, Beck J, Tissler J et al. 2009 Nat. Mater. 8 383–387
  * [4] Ohno K, Joseph Heremans F, Bassett L C, Myers B A, Toyli D M, Bleszynski Jayich A C, Palmstrøm C J and Awschalom D D 2012 Appl. Phys. Lett. 101 082413
  * [5] Beveratos A, Kuhn S, Brouri R, Gacoin T, Poizat J P and Grangier P 2002 Eur. Phys. J. D 18 191–196
  * [6] Aharonovich I, Greentree A D and Prawer S 2011 Nat. Photonics 5 397–405
  * [7] Metsch M H, Senkalla K, Tratzmiller B, Scheuer J, Kern M, Achard J, Tallaire A, Plenio M B, Siyushev P and Jelezko F 2019 Phys. Rev. Lett. 122 190503
  * [8] Rondin L, Tetienne J P, Hingant T, Roch J F, Maletinsky P and Jacques V 2014 Rep. Prog. Phys. 77 056503
  * [9] Jiang L, Hodges J, Maze J, Maurer P, Taylor J, Cory D, Hemmer P, Walsworth R L, Yacoby A, Zibrov A S et al. 2009 Science 326 267–272
  * [10] Fuchs G, Burkard G, Klimov P and Awschalom D 2011 Nat. Phys. 7 789–793
  * [11] Pezzagna S and Meijer J 2021 Appl. Phys. Rev. 8 011308
  * [12] Ruf M, Wan N H, Choi H, Englund D and Hanson R 2021 J. Appl. Phys. 130 070901
  * [13] Casola F, Van Der Sar T and Yacoby A 2018 Nat. Rev. Mater. 3 1–13
  * [14] Schirhagl R, Chang K, Loretz M and Degen C L 2014 Ann. Rev. Phys. Chem. 65 83–105
  * [15] Kehayias P, Levine E, Basso L, Henshaw J, Ziabari M S, Titze M, Haltli R, Okoro J, Tibbetts D, Udoni D et al. 2022 Phys. Rev. Appl. 17 014021
  * [16] Basso L, Kehayias P, Henshaw J, Ziabari M S, Byeon H, Lilly M P, Bussmann E, Campbell D M, Misra S and Mounce A M 2022 Nanotechnology 34 015001
  * [17] Bhaskar M K, Sukachev D D, Sipahigil A, Evans R E, Burek M J, Nguyen C T, Rogers L J, Siyushev P, Metsch M H, Park H et al. 2017 Phys. Rev. Lett. 118 223603
  * [18] Faraon A, Barclay P E, Santori C, Fu K M C and Beausoleil R G 2011 Nature Photonics 5 301–305
  * [19] Riedel D, Söllner I, Shields B J, Starosielec S, Appel P, Neu E, Maletinsky P and Warburton R J 2017 Phys. Rev. X 7 031040
  * [20] Pompili M, Hermans S L, Baier S, Beukers H K, Humphreys P C, Schouten R N, Vermeulen R F, Tiggelman M J, dos Santos Martins L, Dirkse B et al. 2021 Science 372 259–264
  * [21] Yip K Y, Ho K O, Yu K Y, Chen Y, Zhang W, Kasahara S, Mizukami Y, Shibauchi T, Matsuda Y, Goh S K et al. 2019 Science 366 1355–1359
  * [22] Pham L M, Le Sage D, Stanwix P L, Yeung T K, Glenn D, Trifonov A, Cappellaro P, Hemmer P R, Lukin M D, Park H et al. 2011 New J. Phys. 13 045021
  * [23] Guo X, Delegan N, Karsch J C, Li Z, Liu T, Shreiner R, Butcher A, Awschalom D D, Heremans F J and High A A 2021 Nano Letters 21 10392–10399
  * [24] Gallheber B C, Fischer M, Mayr M, Straub J and Schreck M 2018 J. Appl. Phys. 123 225302
  * [25] Challier M, Sonusen S, Barfuss A, Rohner D, Riedel D, Koelbl J, Ganzhorn M, Appel P, Maletinsky P and Neu E 2018 Micromachines 9 148
  * [26] Lee J C, Magyar A P, Bracher D O, Aharonovich I and Hu E L 2013 Diam. Relat. Mater. 33 45–48
  * [27] Piracha A H, Rath P, Ganesan K, Kuhn S, Pernice W H and Prawer S 2016 Nano Lett. 16 3341–3347
  * [28] Aharonovich I, Lee J C, Magyar A P, Buckley B B, Yale C G, Awschalom D D and Hu E L 2012 Adv. Mater. 24 OP54–OP59
  * [29] Piracha A, Ganesan K, Lau D W, Stacey A, McGuinness L P, Tomljenovic-Hanic S and Prawer S 2016 Nanoscale 8 6860–6865
  * [30] Hodges J, Li L, Lu M, Chen E, Trusheim M, Allegri S, Yao X, Gaathon O, Bakhru H and Englund D 2012 New J. Phys. 14 093004
  * [31] Pauls A, Lekavicius I and Wang H 2020 Opt. Express 28 27300–27307
  * [32] Westerhausen M T, Trycz A T, Stewart C, Nonahal M, Regan B, Kianinia M and Aharonovich I 2020 ACS Appl. Mater. Interfaces 12 29700–29705
  * [33] Gaathon O, Hodges J, Chen E, Li L, Bakhru S, Bakhru H, Englund D and Osgood Jr R 2013 Opt. Mater. 35 361–365
  * [34] Evans S 1980 Proc. R. Soc. London, Ser. A 370 107–129
  * [35] Sandstrom R, Ke L, Martin A, Wang Z, Kianinia M, Green B, Gao W b and Aharonovich I 2018 Opt. Commun. 411 182–186
  * [36] Fairchild B A, Olivero P, Rubanov S, Greentree A D, Waldermann F, Taylor R A, Walmsley I, Smith J M, Huntington S, Gibson B C et al. 2008 Adv. Mater. 20 4793–4798
  * [37] Ziegler J F, Ziegler M D and Biersack J P 2010 Nucl. Instrum. Methods Phys. Res. B 268 1818–1823
  * [38] Kehayias P, Henshaw J, Ziabari M S, Titze M, Bielejec E, Lilly M P and Mounce A M 2021 Nucl. Instrum. Methods Phys. Res. B 500 52–56
  * [39] Fairchild B A, Rubanov S, Lau D W, Robinson M, Suarez-Martinez I, Marks N, Greentree A D, McCulloch D and Prawer S 2012 Adv. Mater. 24 2024–2029
  * [40] Sadezky A, Muckenhuber H, Grothe H, Niessner R and Pöschl U 2005 Carbon 43 1731–1742
  * [41] Di Liscia E J, Álvarez F, Burgos E M, Halac E B, Huck H A and Reinoso M E 2013 Mater. Sci. Appl. 04 191–7
  * [42] Orwa J, Nugent K, Jamieson D and Prawer S 2000 Phys. Rev. B 62 5461
  * [43] Jamieson D, Prawer S, Nugent K and Dooley S 1995 Nucl. Instrum. Methods Phys. Res. B 106 641–645
  * [44] Ferrari A C and Robertson J 2004 Philos. Trans. R. Soc. London Ser. A 362 2477–2512
  * [45] Ozawa H, Ishiwata H, Hatano M and Iwasaki T 2018 Physi. Status Solidi 215 1800342
  * [46] Deák P, Aradi B, Kaviani M, Frauenheim T and Gali A 2014 Phys. Rev. B 89 075203
  * [47] Basso L, Bazzanella N, Cazzanelli M and Miotello A 2019 Carbon 153 148–155
  * [48] Bradac C, Gaebel T, Naidoo N, Sellars M, Twamley J, Brown L, Barnard A, Plakhotnik T, Zvyagin A and Rabeau J 2010 Nat. Nanotech. 5 345–349
  * [49] Stacey A, Dontschuk N, Chou J P, Broadway D A, Schenk A K, Sear M J, Tetienne J P, Hoffman A, Prawer S, Pakes C I et al. 2019 Adv. Mater. Interfaces 6 1801449
  * [50] Magyar A P, Lee J C, Limarga A M, Aharonovich I, Rol F, Clarke D R, Huang M and Hu E L 2011 Appl. Phys. Lett. 99 081913
  * [51] Jin K, Velisa G, Xue H, Yang T, Bei H, Weber W J, Wang L and Zhang Y 2019 J. Nucl. Mater. 517 9–16


# On a conjecture of Erdős

Yong-Gao Chen and Yuchen Ding  (Yong-Gao Chen) School of Mathematical Science,
Nanjing Normal University, Nanjing 210023, People’s Republic of China
<EMAIL_ADDRESS>(Yuchen Ding) School of Mathematical Science, Yangzhou
University, Yangzhou 225002, People’s Republic of China<EMAIL_ADDRESS>
###### Abstract.

Let $\mathcal{P}$ denote the set of all primes. In 1950, P. Erdős conjectured
that if $c$ is an arbitrarily given constant, $x$ is sufficiently large and
$a_{1},\dots,a_{t}$ are positive integers with

$a_{1}<a_{2}<\cdot\cdot\cdot<a_{t}\leqslant x,~{}t>\log x,$

then there exists an integer $n$ so that the number of solutions of
$n=p+a_{i}$ $(p\in\mathcal{P},1\leq i\leq t)$ is greater than $c$. In this
note, we confirm this old conjecture of Erdős.

###### Key words and phrases:

Representation function; Primes

###### 2010 Mathematics Subject Classification:

Primary 11A41; Secondary 11A67.

*Corresponding author

##  1\. Introduction

Let $\mathcal{P}$ denote the set of all primes. In 1950, Erdős [5] made the
following anecdotal conjecture:

Erdős Conjecture. Let $c$ be any constant and $x$ sufficiently large,

$a_{1}<a_{2}<\cdot\cdot\cdot<a_{t}\leqslant x,~{}t>\log x.$

Then there exists an integer $n$ so that the number of solutions of
$n=p+a_{i}$ $(p\in\mathcal{P},1\leq i\leq t)$ is greater than $c$.

Erdős [5] himself proved this conjecture for the case $a_{i}=2^{i}$, which
gives an affirmative answer to a question of Turán. In a former note [3], the
second author proved this conjecture for the case $a_{i}\mid a_{i+1}$ with its
quantitative form, which is a slight generalization of Erdős’ result. In a
subsequent note, the second author and Zhou [4] proved the conjecture for the
case $a_{i}=2^{p_{i}}$, where $p_{i}$ is the $i$-th prime. This case was
conjectured by the first author [2] years ago. Shortly after, the authors of
the present note recognized that the complete proof of Erdős’ conjecture
actually follows directly from a new achievement of the distributions of the
primes established by Maynard–Tao [7, 9]. We keep record here as the closure
of this longstanding conjecture.

In this note, the following general results are proved. The Erdős Conjecture
follows from Corollary 1.2.

######  Theorem 1.1.

For any $\ell$ distinct integers $a_{1},\dots,a_{\ell}$, there are infinitely
many positive integers $n$ such that the number of solutions of $n=p+a_{i}$
$(p\in\mathcal{P},1\leq i\leq\ell)$ is greater than

$\frac{1}{8}\log\ell-1.6.$

From Theorem 1.1, we immediately have the following corollaries:

######  Corollary 1.2.

Let $x\geq 2$ and

$a_{1}<a_{2}<\cdot\cdot\cdot<a_{t}\leqslant x,~{}t>\log x.$

Then there exists infinitely many integers $n$ so that the number of solutions
of $n=p+a_{i}$ $(p\in\mathcal{P},1\leq i\leq t)$ is greater than

$\frac{1}{8}\log\log x-1.6.$

######  Corollary 1.3.

Let $\mathcal{A}=\\{a_{i}\\}_{i=1}^{\infty}$ be an infinite set of integers
and let

$f_{\mathcal{A}}(n)=\\#\\{(p,a):n=p+a,p\in\mathcal{P},a\in\mathcal{A}\\},$

then

$\limsup_{n\to+\infty}f_{\mathcal{A}}(n)=+\infty.$

##  2\. Proofs

A set $\\{b_{1},\dots,b_{k}\\}$ is called an admissible set if there is no a
fixed integer $d>1$ such that $d\mid(n+b_{1})\cdots(n+b_{k})$ for all integers
$n$. It is equivalent that $\\{b_{1},\dots,b_{k}\\}$ modulo $p$ occupies at
most $p-1$ residues. We begin with the following deep result for the
distribution of the primes due to Maynard–Tao.

######  Lemma 2.4.

[6, Theorem 6.2] For any given integer $m\geqslant 2$, let $k$ be a positive
integer with $k\log k>e^{8m+4}$. For any admissible set
$\\{b_{1},\dots,b_{k}\\}$, there are infinitely many integers $n$ such that at
least $m$ of $n+b_{1},\dots,n+b_{k}$ are prime numbers.

######  Lemma 2.5.

[1, Lemma 3] We have

$\prod_{3\leq p\leq x}\left(1-\frac{1}{p}\right)^{-1}\leq 0.923\log x,\quad
x\geq 74,$

where the product is taken over all primes $p$ with $3\leq p\leq x$.

###### Proof of Theorem 1.1.

If $\ell\leq e^{12}$, then

$\frac{1}{8}\log\ell-1.6\leq 0,$

and Theorem 1.1 is trivial. In the following, we assume that $\ell>e^{12}$.

Let $p_{i}$ be the $i$-th prime. Assume that $a_{1},\dots,a_{\ell}$ are $\ell$
distinct integers. For $p_{1}$, one of residues modulo $p_{1}$ contains at
most $\lfloor\ell/p_{1}\rfloor$ of $a_{1},\dots,a_{\ell}$. So at least
$\ell-\lfloor\ell/p_{1}\rfloor$ of $a_{1},\dots,a_{\ell}$ occupy at most
$p_{1}-1$ residues modulo $p_{1}$. Let $\ell_{0}=\ell$ and
$\ell_{1}=\ell-\lfloor\ell/p_{1}\rfloor$. Without loss of generality, we
assume that $a_{1},\dots,a_{\ell_{1}}$ occupy at most $p_{1}-1$ residues
modulo $p_{1}$. Similarly, without loss of generality, we may assume that
$a_{1},\dots,a_{\ell_{2}}$ occupy at most $p_{2}-1$ residues modulo $p_{2}$,
where $\ell_{2}=\ell_{1}-\lfloor\ell_{1}/p_{2}\rfloor$. Continuing this
process, at the $t$-th step, we may assume that $a_{1},\dots,a_{\ell_{t}}$
occupy at most $p_{t}-1$ residues modulo $p_{t}$, where
$\ell_{t}=\ell_{t-1}-\lfloor\ell_{t-1}/p_{t}\rfloor$. Since
$\ell\geq\ell_{1}\geq\cdots$ and $p_{1}<p_{2}<\cdots$, there exists $t$ with
$\ell_{t}<p_{t+1}$. Let $s$ be the least integer with $\ell_{s}<p_{s+1}$. It
is clear that $\\{a_{1},\dots,a_{\ell_{s}}\\}$ is an admissible set. Let $m$
be largest integer with $\ell_{s}\log\ell_{s}>e^{8m+4}$. If $m\geq 2$, then by
Lemma 2.4, there are infinitely many integers $n$ such that at least $m$ of
$n-a_{1},\dots,n-a_{\ell_{s}}$ are prime numbers. Since there are infinitely
many primes, it follows that there are infinitely many integers $n$ such that
at least one of $n-a_{1},\dots,n-a_{\ell_{s}}$ is prime number. So the
conclusion is also true for $m\leq 1$.

Now we establish an explicit relation between $\ell$ and $m$.

Since

$\ell_{i+1}=\ell_{i}-\left\lfloor\frac{\ell_{i}}{p_{i+1}}\right\rfloor\geq\ell_{i}-\frac{\ell_{i}}{p_{i+1}}=\ell_{i}\left(1-\frac{1}{p_{i+1}}\right),\quad
i=0,1,\dots,$

it follows from the definition of $s$ that

$\displaystyle p_{s+1}$ $\displaystyle>$
$\displaystyle\ell_{s}\geq\ell_{s-1}\left(1-\frac{1}{p_{s}}\right)\geq\cdots$
(2.1) $\displaystyle\geq$
$\displaystyle\ell\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{s}}\right)$
$\displaystyle>$ $\displaystyle
e^{12}\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{s}}\right).$

Noting that

$p_{101}<e^{12}\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{100}}\right),$

for $1\leq i\leq 100$ we have

$\displaystyle p_{i+1}\leq p_{101}$ $\displaystyle<$ $\displaystyle
e^{12}\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{100}}\right)$
$\displaystyle\leq$ $\displaystyle
e^{12}\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{i}}\right).$

It follows from (2.1) that $s>100$. So $p_{s}\geq p_{101}=547$. Thus, by Lemma
2.5,

$\ell_{s}\geq\ell\left(1-\frac{1}{p_{1}}\right)\cdots\left(1-\frac{1}{p_{s}}\right)\geq\frac{\ell}{2}\cdot\frac{1}{0.923\log
p_{s}}=\frac{\ell}{1.846\log p_{s}}.$

By the definition of $s$, $p_{s}\leq\ell_{s-1}$. Thus,

$\ell_{s}\geq\left(1-\frac{1}{p_{s}}\right)\ell_{s-1}\geq\left(1-\frac{1}{p_{s}}\right)p_{s}=p_{s}-1.$

It follows that

$\ell_{s}\geq\frac{\ell}{1.846\log p_{s}}\geq\frac{\log 546}{1.846\log
547}\frac{\ell}{\log(p_{s}-1)}>\frac{0.54\ell}{\log\ell_{s}}.$

So $\ell_{s}\log\ell_{s}\geq 0.54\ell$. In view of the definition of $m$,

$e^{8m+12}\geq\ell_{s}\log\ell_{s}\geq 0.54\ell.$

So

$m\geq\frac{1}{8}\log\ell-\frac{12}{8}+\frac{\log
0.54}{8}>\frac{1}{8}\log\ell-1.6.$

This completes the proof of Theorem 1.1. ∎

###### Proof of Corollary 1.2.

Assume that

$1\leq a_{1}<\cdots<a_{t}\leq x,\quad t>\log x.$

By Theorem 1.1, there are infinitely many positive integers $n$ such that the
number of solutions of $n=p+a_{i}$ $(p\in\mathcal{P},1\leq i\leq t)$ is
greater than

$\frac{1}{8}\log t-1.6>\frac{1}{8}\log\log x-1.6.$

This completes the proof of Corollary 1.2. ∎

###### Proof of Corollary 1.3.

By Theorem 1.1, there is a positive integers $n$ such that the number of
solutions of $n=p+a_{i}$ $(p\in\mathcal{P},1\leq i\leq\ell)$ is greater than
$\frac{1}{8}\log\ell-1.6$. That is,
$f_{\mathcal{A}}(n)\geq\frac{1}{8}\log\ell-1.6$. Now Corollary 1.3 follows
immediately. ∎

##  3\. Remarks

It is known that there is a positive proportion of positive odd numbers that
can be represented as $p+2^{k}$ with $k\in\mathbb{N}$ and $p\in\mathcal{P}$
(See Romanoff [8] ) and there is an arithmetical progression of positive odd
numbers none of which can be represented as $p+2^{k}$ with $k\in\mathbb{N}$
and $p\in\mathcal{P}$ (see Erdős [5]).

Since one can take $k\leq cm^{2}e^{4m}$ for some positive constant $c$ in
Lemma 2.4 (see [6] ), it follows that $\frac{1}{8}$ in Theorem 1.1 and
Corollary 1.2 can be improved to any constant less than $\frac{1}{4}$ for all
sufficintly large $\ell$ and $x$.

## Acknowledgments

The first author was supported by the National Natural Science Foundation of
China, Grant No. 12171243. The second author was supported by the Natural
Science Foundation of Jiangsu Province of China, Grant No. BK20210784. He was
also supported by the foundations of the projects “Jiangsu Provincial
Double–Innovation Doctor Program”, Grant No. JSSCBS20211023 and “Golden Phenix
of the Green City–Yang Zhou” to excellent PhD, Grant No. YZLYJF2020PHD051.

## References

  * [1] Y.-G. Chen, X.-G. Sun, On Romanoff’s constant, J. Number Theory, 106 (2004), 275–284.
  * [2] Y.-G. Chen, Romanoff theorem in a sparse set, Sci. China Math., 53 (2010), 2195–2202.
  * [3] Y. Ding, A note on a conjecture of Erdős, (preprint).
  * [4] Y. Ding, G.-L. Zhou, Some application of the admissible sets, (preprint).
  * [5] P. Erdős, On the integers of the form $2^{k}+p$ and some related problems, Summa Brasil. Math., 2 (1950), 113–123.
  * [6] A. Granville, Primes in intervals of bounded length, Bull. Amer. Math. Soc., 52 (2015), 171–222.
  * [7] J. Maynard, Small gaps between primes, Ann. of Math., 181 (2015), 383–413.
  * [8] N.P. Romanoff, Über einige Sätze der additiven Zahlentheorie, Math. Ann. 57 (1934) 668–678.
  * [9] D.H.J. Polymath, Variants of the Selberg sieve, and bounded intervals containing many primes, Res. Math. Sci., 1 (2014), 1–83.