\section{Introduction}
Ultra-wideband optical transmission systems allow increased fiber link capacities, and are being used by network operators aiming to fully exploit their existing single-mode fiber (SMF) infrastructure. Currently, most optical communication systems are operating over the C-band (5~THz wide), with some high speed links being operated over the entire C+L band (10~THz wide). Experimental demonstrations of C+L band systems transmitting over transoceanic distances date back to 2013 \cite{Salsi2013transmission}. Ever since, the capacity of transoceanic C+L systems gradually increased as shown in Fig. 1 (after \cite[Fig. 1]{Ionescu2019Transmission}), mainly driven by advanced coded modulation and digital signal processing. Modern C+L band systems achieve an approximate twofold increase compared to their C-band only counterparts. While C+L band systems are a commercial reality and the demand for fiber capacity is still growing, expanding transmission bandwidths to the entire S+C+L-band (20~THz wide) is increasingly explored. Enabled by ultra-wideband amplification technologies, such as thulium-doped fiber amplifiers (TDFA) and semiconductor optical amplifiers (SOA), first S+C+L band transmission systems were recently demonstrated over a continuous bandwidth of 12.5~THz \cite{Renaudier_2017_f1c} and discontineous bandwidths of 16.4~THz \cite{Hamaoka_2018_usf} and 15.1~THz \cite{Puttnam_2020_059}. However, to predict the capacity benefits of the S-band in mature and commercially available S+C+L band systems in the future, analytical models can be utilized. Recently, the formalism of the Gaussian Noise (GN) model was extended to include inter-channel stimulated Raman scattering (ISRS), enabling accurate predictions of ultra-wideband transmission performance \cite{Roberts_17_cpo, Cantono_2018_oti, Semrau_2018_tgn}. Furthermore, a modulation format dependent closed-form formula was introduced enabling large scale system optimization and capacity estimations for arbitrary modulation formats~\cite{Semrau_2019_aca, Semrau_2019_amf,Semrau_2019_mfd}.
\begin{wrapfigure}[16]{r}{0.515\textwidth}
   \centering
    \includegraphics[]{figure0.pdf}
\vspace{-0.1cm}
\caption{Records of transoceanic transmission experiments over $>5900$~km. Also shown, modelling predictions as a result of this work using the ISRS GN model.}
\label{fig:capacity}
\end{wrapfigure}
\
In this paper, the throughput gains of extending the transmission window to the entire S+C+L band are analytically estimated using the modulation format dependent ISRS GN model. Two scenarios are studied; 1) using ideal dynamic gain equalization at each amplification stage, and 2) considering non-ideal gain equalization. The latter leads to ISRS accumulation across multiple fiber spans. Two distinct, geometrically shaped modulation formats and a finite set of available code rates across the optical spectrum are considered.
\section{Analytical performance prediction}
To estimate the throughput of the system under test, we use a modulation format dependent, closed-form approximation of the ISRS GN model \cite{Semrau_2019_aca, Semrau_2019_amf}. The formula was derived assuming a triangular approximation of the Raman gain spectrum, valid up to 15~THz optical bandwidth. However in \cite{Semrau_2019_mfd}, we proposed a semi-analytical approach that enables the formula to be applied to bandwidths beyond 15~THz. The approach consists of matching the actual signal power profile to its first-order approximation (whose parameter are then used in the formula). The methodology yields very fast, but yet accurate, performance estimations and enables optimization problems with large sample space, such as launch power optimization.
We consider 364$\times$50~GBd channels, centered at 1540~nm and transmitted over a transoceanic distance of 7000~km (100$\times$70~km). A transoceanic distance was chosen to ease the comparison to experimental demonstrations. Spectral gaps of 10~nm and 5~nm were considered between the respective transmission bands. A Corning\textsuperscript{\textregistered} SMF-28\textsuperscript{\textregistered} ULL fiber was assumed with dispersion and nonlinearity parameter as $D=18\frac{\text{ps}}{\text{nm}\cdot \text{km}}$, $S=0.067\frac{\text{ps}}{\text{nm}^2\cdot\text{km}}$ and $\gamma=1.2\frac{\text{1}}{\text{W}\cdot\text{km}}$ (see \cite{Semrau_2019_mfd} for more details). The noise figures of the EDFA/TDFA were 7~dB, 4~dB, and 6~dB in the S, C and L-band, respectively. The spectral launch power distribution was optimized using a particle swarm algorithm with subsequent gradient descent. Two distinct, geometrically shaped modulation formats were considered. In the S-band a 16-QAM format was transmitted (optimized for 7~dB SNR with excess kurtosis $\Phi=-0.49$), whereas in the C+L band a 64-QAM format was transmitted (optimized for 11~dB SNR with excess kurtosis $\Phi=-0.32$). Once the signal-to-noise ratio (SNR) was calculated, a finite set of code rates was selected. The code rates were optimized to maximize the total throughput and were upper bounded by the normalized generalized mutual information (GMI) of the respective channels.
Two scenarios were considered: 1) ideal gain equalization, where the transmitted launch power was recovered at each span, and 2) non-ideal gain equalization, where after each span only 50\% of the ISRS power transfer is compensated and after every fifth span the transmitted launch power is ideally recovered. In the latter case, the ISRS power transfer is accumulating across multiple spans. The optimized SNR and launch power distribution are shown in Fig. 2a) and b), respectively. The inset in Fig. 2a) shows the considered modulation formats. The GMI and the total throughput depending on the number of code rates for both cases are shown in Fig. 2c). For comparison, the analysis was repeated for a C+L band system with the achieved throughput shown in Fig. 1.
\begin{figure*}[t!]
   \centering
    \includegraphics[]{figure1.pdf}
\vspace{-1ex}
\caption{\small The a) SNR and b) launch power after $\text{100}\times\text{70~km}$ transmission for ideal and non-ideal gain equalization. The throughput and GMI vs. the number of code rates are shown in c). The insets in a) show the used modulation formats in the respective transmission bands.}
\vspace{-3ex}
\label{fig:NLI}
\end{figure*}
The results in this paper show that current throughputs can be increased from 74.59~Tbit/s to around 119.5~Tbit/s (achieving 119.2~Tbit/s using 6 code rates) when using the S-band for fiber transmission. This can be achieved by launch power optimization, using two distinct modulation formats for the S and C+L band and only a few code rates. It was also shown that non-ideal gain equalization reduces the throughput to around 112.3~Tbit/s, emphazing the importance of ISRS mitigation. Although the analysis in this paper was carried out for transoceanic distances, the conclusions carry over to terrestial distances.

\section{Conclusions}
The benefits of using the S-band in optical communication systems were studied, indicating throughput gains of around 60\% compared with C+L band transmission. The capacity gain can be unlocked using different modulation formats in the respective transmission bands, spectrally optimized launch powers and dynamic gain equalization.

\bibliographystyle{IEEEtran}


\subsection*{Introduction}
	Topological phase, possessing intriguing bulk and edge properties, plays a fundamental role in understanding matter~\cite{Hall1,review1,review2,review3,review4} and constructing artificial devices. It displays an extraordinary robustness to smooth changes in material parameters or disorders and endows the system with inherent protection~\cite{review1, review2}. In addition to fundamental physics, the extraordinary robustness of topological phase has found promising applications for inherently fault-tolerant quantum simulation and quantum computing~\cite{computing}. In the past decades, topological robust phases have been widely investigated in various systems including condensed-matter~\cite{review1, review2}, ultracold matter~\cite{Cold1, Cold2, Cold3}, phononic~\cite{Topo_weyl_1,XBY4,XBY1} and photonic systems~\cite{review3, review4, photon1, photon2, photon3}.

	Compared with other information carriers, photon is easier to be created and detected, and the topological photonics is rapidly developed in recent years. Together with the three-dimensional femtosecond laser direct writing technique~\cite{fabri_1}, we are able to construct photonic waveguide lattice on a photonic chip in a physically scalable and structurally designable fashions~\cite{LargestChip, 3Dscalable, PIT_Gap}. The constructed photonic waveguide array provides a convenient platform for investigating the topological phase, especially, for the higher-order topological insulators (HOTIs) and topological crystalline insulators (TCIs) in a high complexity~\cite{TCI}.

	HOTIs have been recently proposed as a novel topological phase of matter with unconventional bulk-boundary correspondence~\cite{HOTI1, HOTI2, HOTI3, HOTI4, HOTI5, HOTI7, HOTI8, HOTI9, HOTI10, HOTI12, HOTI13, HOTI14, HOTI15, corner, XBY3, QP1, quadrupole, C3}. An $N$th-order topological insulator has topological boundary state with codimension $N$~\cite{CO1,CO2}. Generally, there are two ways to achieve HOTIs. The first one is to use topological multipole insulators which were firstly proposed by W. Benalcazar et al.~\cite{HOTI1,QP1} and lately experimentally realized in mechanics~\cite{HOTI2}, microwaves~\cite{HOTI3}, electric circuits~\cite{HOTI4}, phononics~\cite{HOTI15} and coupled-resonator optical waveguides quadrupole~\cite{quadrupole}. The other approach is to utilize topological crystalline insulators with quantized bulk polarization which was theoretically proposed by considering tight-binding models~\cite{HOTI8,CO1} and experimentally achieved in photonics~\cite{HOTI9, HOTI10, XBY3, C3} and phononics~\cite{HOTI12, HOTI13, HOTI14}. However, the investigations of the intriguing lower dimensional topological boundary states in HOTIs down to single-particle level as well as their quantum characters have not yet been explored. The impact of topology on quantum entanglement for entangled photon states also attracts the interest, which inspire a fascinating and elegant combination of topological photonics and quantum photonics.

	Here, we theoretically and experimentally investigate the higher-order topological phases and the topological protection to quantum superposition and entanglement in two-dimensional photonic lattice via single-photon dynamics. In contrast to all-dielectric photonic crystals, the zero-dimensional corner state generated in the second-order topological phases is more sensitive to the crystalline symmetry of photonic lattice. We experimentally demonstrate two ways of identifying the crystal-symmetry-dependent topological corner states in the photonic lattice, with single-site and superposition-state injection. The topological corner state is found being able to preserve single-photon superposition state and the quantum entanglement against diffusion-induced decoherence, representing an intrinsic topological protection mechanism for multi-partite quantum states.\\

	\subsection*{Experimental implementation and results}
	\paragraph*{Integrated topological lattice.} In our experiment, we fabricate the higher-order topological insulators in alkaline earth boro-aluminosilicate glass using femtosecond laser direct writing technique (see Methods)~\cite{fabri_1}. We integrate various samples in a photonic chip, as shown in Fig.~\ref{f1}(a). The constructed two-dimensional lattice contains $8\times8$ sites and the evolution distances vary from 10 to 30 mm with step of 5 mm. Our photonic lattice mimics the two dimensional (2D) extended Su-Schrieffer-Heeger (SSH) model~\cite{SSH1}, the corresponding Hamiltonian in the momentum space can be expressed as
	\begin{equation}
		\mathcal{H}(\bm{k})=
		\left(
		\begin{array}{cccc}
			0 & h_{12} & h_{13} & 0\\
			h_{12}^\ast & 0 & 0 & h_{24}\\
			h_{13}^\ast & 0 & 0 & h_{34}\\
			0 & h_{24}^\ast & h_{34}^\ast & 0 \\
		\end{array}
		\right),
		\label{eq:1}
	\end{equation}
	where $h_{12}=t^x_a+t^x_b\mathrm{exp}(\mathrm{i}k_x)$, $h_{13}=t^y_a+t^y_b\mathrm{exp}(\mathrm{-i}k_y)$, $h_{24}=t^y_a+t^y_b\mathrm{exp}(-\mathrm{i}k_y)$, $h_{34}=t^x_a+t^x_b\mathrm{exp}(\mathrm{i}k_x)$, and ${\bm k}=(k_x, k_y)$. The parameter $t^i_a$ ($t^i_b$) with $i=x,y$ represents the intra-cell (inter-cell) coupling strength along $i$-direction. The coupling strength is modulated by the separation between two sites, as shown in Fig.~\ref{f1}(b), the corresponding separation distance for $t^i_a$ ($t^i_b$) is $d^i_a$ ($d^i_b$). For our waveguide array, the coupling between nearest-neighbour sites are positive which means $t^i_j>0$ for $i=x,y$ and $j=a,b$ and the coupling between next-nearest-neighbour (or higher-order-neighbour) sites are exponentially suppressed~\cite{2ndC} which ensures the validity of the tight-binding approximation.

	By diagonalizing the Hamiltonian in Eq.(\ref{eq:1}), we can obtain the band structure of the photonic lattice [see Figs.\ref{f1}(c)-\ref{f1}(h)]. The competition between $t^i_a$ and $t^i_b$ with $i=x, y$ determines the existence of the band gap where the band gap closes at the topological phase transition point $t^i_a=t^i_b$. On the other hand, if $t^x_j \neq t^y_j$ ($t^x_j = t^y_j$) for $i=a,b$, the system described by the above Hamiltonian always have $C_2$ ($C_4$) rotation symmetry. The band structures of $C_4$ and $C_2$ symmetric lattices with different coupling parameters are presented in Figs.\ref{f1}(c)-\ref{f1}(e) and Figs.\ref{f1}(f)-\ref{f1}(h) respectively, where the corresponding relationship are $t^i_a>t^i_b$, $t^i_a=t^i_b$, and $t^i_a<t^i_b$. A band inversion process has been observed in both cases, indicating a topological phase transition when $t^i_a=t^i_b$.

	\paragraph*{Topological corner index and filling anomaly.} The topological photonic lattices with $C_4$ and $C_2$ crystalline symmetry both are in the second-order topological phases with 1D topological edge states and 0D topological corner states (see Methods). To investigate the corner physics of the second-order topological insulator (SOTI), we explore the rotation symmetry group representations at high symmetry points in Brillouin zone~\cite{TCI1,HOTIC}. Due to the symmetries, the eigenstates of the Hamiltonian can be chosen as the common eigenstates of the rotation operators $\hat{R}_n$ ($n=2,4$) with corresponding eigenvalues $\Pi^{(n)}_p=e^{2\pi i(p-1)/n}$ with $p=1,2,...,n$. By comparing the rotation eigenvalues at high symmetry points to those of atomic insulators, we can determine whether the system is topologically non-trivial. Specifically, we can define the topological index as
	\begin{equation}
		[\Pi^{(n)}_p]=\#\Pi^{(n)}_p-\#\Gamma^{(n)}_p,
		\label{index}
	\end{equation}
	where $\#\Pi^{(n)}_p$ is the number of bands below the band gap with rotation eigenvalues $\Pi^{(n)}_p$. $\Pi^{(n)}$ stand for high symmetric point $X$, $M$ and $\Gamma$ in $C_n$ symmetric systems. When $[\Pi^{(n)}_p]$ is non-zero, the system is a topological insulator (denoted as the obstructed atomic insulator).

	Theoretically the indices defined in Eq.\ref{index} can fully characterize the higher-order topological properties of the systems. However, considering the time-reversal symmetry and the fact that the number of bands below the gap is constant across the Brillouin zone, these indices are not independent to each other, thus we can drop the redundant indices. For $C_4$ symmetric lattices, the indices are $[X_1]$, $[M^{(4)}_1]$ and $[M^{(4)}_2]$ while for $C_2$ symmetric lattices, the indices are $[X_1]$, $[Y_1]$ and $[M^{(1)}_1]$. This approach is based on the theory of topological crystalline insulators~\cite{HOTIC} and similar to the recently proposed topological symmetry indicators~\cite{SI1,SI2,SI3}.

	We further determine the values of these topological indices by investigating the coupling configurations of lattice. We find that it is $\frac{1}{4}$ fractionalized corner states at each of four corners for our $C_4$ symmetric lattice, but $\frac{1}{2}$ fractionalized corner states at each of four corners for $C_2$ symmetric lattice (see Methods for details). The topological corner index captures the corner physics more precisely and directly than the bulk polarization. Moreover, it provides the fractionalization of the photonic eigenstates. Due to the Abelian additive structure of TCI, we can theoretically construct TCIs with other topological corner index by using these primitive generators and even for the fragile TCIs that do not admit Wannier representations \cite{HOTIC}.

	\paragraph*{Chiral symmetry and the separation between corner states and bulk states.} Different from the all-dielectric photonic crystals~\cite{XBY3}, the energy levels between corner states and bulk states of our system is able to manifest the difference between $C_2$ and $C_4$ rotation symmetry. In 1D SSH model~\cite{SSH_1,SSH_2,SSH_3}, due to the staggered coupling strengths between inter-cell and intra-cell sites, there is a sublattice symmetry which is also the chiral symmetry of the Hamiltonian since it satisfies $\Gamma H(\bm{k})\Gamma^{-1}=-H(\bm{k})$, where $\Gamma$ presents operation of chiral symmetry. It restricts the band structures (or more precisely, the eigenvalues) to be symmetric with respect to ``zero energy". In 2D cases, due to the exponentially suppressing of the higher-order couplings, the chiral symmetry is preserved in $C_4$ symmetric lattice while being broken in $C_2$ symmetric lattice. Therefore, in $C_4$ symmetric case, the corner states are pinned on the zero-energy level and embedded into the bulk states. In fact, these corner states are actually the higher-order version of bound states in the continuum, the details analysis can be found in Supplementary Materials. On contrary, in $C_2$ symmetric case, the corner states are not restricted to the zero-energy level and can be separated from the bulk states and edge states. The above analysis can be confirmed by our following experimental observation.

	\paragraph*{Corner states in finite lattice.} For the finite lattice, the edge states and the corner states manifest due to the dimensional hierarchy of topological phases. We show the spectrum of the $C_4$ and $C_2$ symmetric lattices in Figs.\ref{f2}(a)-\ref{f2}(d) respectively. For the finite $C_4$ symmetric lattice, we show the spectrum as the function of $t_a/t_b$ varying from -1 to 1 in Fig.~\ref{f2}(a), and the detailed spectrum for the case of $t_a/t_b=0.1$ in Fig.~\ref{f2}(b). We further show the spatial distribution of the zero-energy corner modes and the edge modes in the insets (i)-(iv) in Fig.~\ref{f2}, which can be defined by
	\begin{equation}
		D_n(E)=\sum_{m}\delta(E-E_m)|\varphi_n^{(m)}|^2,
	\end{equation}
	where $E_m$ is the energy of the $m$th eigenstate $\varphi_n^{(m)}$. As the spatial distribution shown, the photon is confined at the four or two corners of lattice under the norm of the zero-energy corner states [see insets (i-ii)]. The edge states imply that the photon occupies the boundaries of the lattice with high probability [see insets (iii)]. In contrast, the photon only distributes in the bulk of the lattice for the bulk states [see insets (iv)].

	For the finite $C_2$ symmetric lattice, the spectrum behaves differently with the $C_4$ symmetric lattice. There are no zero-energy bulk states near the corner states, and no eigen-diagonal corner state among the four corner eigenstates as the $C_2$ symmetric lattice. Moreover, for the $C_2$ lattice, the corner states are two-fold degenerate as shown in insert of Fig.~\ref{f2}(b), which is characterized by $Q^{(2)}_c=\frac{1}{2}$. The edge states have the same symmetry as the lattice, as a result, there are two kind eigen edge states in $C_2$ symmetric lattice [see insets (vi-vii)].

	\paragraph*{Identifying the corner states.} As analyzed above, the corners state behaves as that the single-photon probability distribution mainly occupies on the corners of lattice. According to the principle of quantum superposition, if we inject the single photon into the lattice from one of the corners, only the corner states mainly maintain the evolution of single photon in lattice, though all the states are excited. In this case, the single photon will be confined in the excited corner of lattice, which is the result of the superposition of all the zero-energy corner states. Such evolution is stable due to that these corner states are robust to the structure disorder and are decoupled from the bulk bands.

	In experiment, the parameters are adopted as $d^x_a=d^y_a=22\ \mu$m, $d^x_b=d^y_b=9\ \mu$m for the $C_4$ symmetric lattice. We inject the single photons into each corner and capture the single-photon distribution probability after different evolution distance varying from 10 to 30 mm with a step of 5 mm [see Fig.~\ref{f2}(e) and Methods for details]. The measured single-photon distribution probability with evolution distance of 10 mm are shown in Fig.~\ref{f2}(f), where the single photon almost only occupies on the excited site. To quantify the localization of the outgoing single-photons distribution, we define the generalized return probability as $\xi=\sum_{k-w}^{k+w}I_i/\sum_{1}^{n}I_i$, where the $n$ is the site number of lattice. The quantity $\xi$ quantifies the probability of the single photon remaining within a small width $w$ from the injected site $k$. For the corner states, the photon is expected to be localized in only the corner site, then the accumulated width $w=0$. Such that the return probability of corner state is $\xi=\sum_{k}I_i/\sum_{1}^{n}I_i$. In our experiment, as shown in Fig.~\ref{f2}(g), the return probability of corner state $\xi$ maintains in high value close to 1 and doesn't change with the increase of evolution distance. We modulate the $d^x_b$ varying from 11 to 14 $\mu$m, the corresponding $t_a/t_b$ varies from 0.08 to 0.22, and keep other parameters unchanged. Five $C_2$ symmetric lattices are integrated in one photonic chip. The measured results are shown in Fig.~\ref{f2}(f)-\ref{f2}(g), which accord to the $C_4$ symmetric lattice as analyzed in theory.

	\paragraph*{Compatibility of corner state and single-photon superposition state.} As we known, the quantum superposition, as the most fundamental principle of quantum mechanics, allows the single photon be in more than one state simultaneously, which is called single-photon superposition state. For example, if there are four configurations labeled by $\ket{\psi_i}$, $i=1,2,3,4$, the most general single-photon superposition state would be

	\begin{equation}
		\ket{\psi}=\sum_{1}^{4}c_i\ket{\psi_i},
	\end{equation}
	where the coefficients are complex numbers describing the probability of each configuration. With the help of such a single-photon superposition state, we are able to determinately identify the corner states without exciting all states.

	As shown in Fig.~\ref{f3}(a)-\ref{f3}(b), we obtain the single-photon superposition state as $\ket{\psi}=\frac{1}{2}\sum_{1}^{4}\ket{\psi_i}$ by injecting a single photon into a 3D 1$\times$4 photonic coupler [see Methods for details]. Subsequently, the prepared single photon is split into four corners of lattice to generate the desired superposition states, and the distribution probability of photon in lattice is identical to the zero-energy corner state. In the language of quantum mechanics, the system now is in the eigenstate and the single-photon superposition state in lattice will be maintained due to the orthogonality among eigenstates. With the same parameters, five kinds lattice are fabricated and integrated with the 3D 1$\times$4 photonic coupler on one chip. In Fig.~\ref{f3}(c)-\ref{f3}(h), we show the measured single-photon distribution probability and the return probability. The output probability distribution of single photon follows the distribution of corner states and is robust under the change of structure parameters of the systems.

	 Meanwhile, the result also implies that the evolution of the single-photon superposition state is also protected by the topological corner mode against the structure disorder and diffusion induced photon loss. The successful observation of topological protected single-photon superposition state here adds one more key element into the toolbox of quantum topological photonics. Together with the previous work demonstrating the topologically protected quantum states and source~\cite{Qutop_1photon,Qutop_source,Qutop_interfere,Qutop_2photon_1,Qutop_2photon_2,Qutop_entanglement}, the robust features against photon loss induced by diffusion and the structure disorder promise that the topological protected single-photon superposition state is able to provide a valuable resource for quantum information processing.

	\paragraph*{Protecting the entanglement.} It is interesting to examine the impact of topological phase on quantum entanglement for entangled photon states~\cite{Qutop_2photon_2, Qutop_entanglement}. In our experiment, we inject the entangled photons into the corner state of the photonic lattice and reconstruct the density matrix of the output states. We calculate the concurrence and purity of the measured two-qubit entangled state to quantify the entanglement of the output photons, where the high concurrence and purity indicate the high entanglement quality and quantum state respectively. As shown in Fig.~\ref{F_Entang}, the entangled photons in the corner state preserve the high concurrence and purity beyond 90\%, even after introducing the disorder. In contrast, the entanglement in the trivial cases tends to corruption, which concurrence and purity drop to lower than 90\%. Especially, we fabricate a large-scale two-dimension uniform lattice with $21\times21$ sites to value the entanglement in the quantum walk. The concurrence and purity of the entangled state in the quantum walk drop to lower than 80\% after evolution distance of $z = 11$ mm, which lose the possibility of in practical applications.

	\paragraph*{The finite gap effect.} As the spectrum shown in Fig.~\ref{f4}(a), the degenerate corner modes will transit to non-degenerate bulk states with the increase of $t_a/t_b$ for the topological phase. This is the consequence of the finite gap effect. Under the same size of lattice, when $t_a/t_b$ approaches to $1$, the gap size between the first band and second band (or the third band and the forth band) approaches to zero. Since the gap size is proportional to the mass term in the effective Hamiltonian which is exponentially proportional to the decay rate of the interface states away from the corners, the photon in corner states will tend to occupy the boundaries and the bulk [see Fig.~\ref{f4}(b)]. Due to the overlap of the evanescent modes between nearest-neighbour waveguides, the photon distribution will evolute to all corners. A 1/4 fractional probability of photon on the four corners can be observed even when we inject a photon in a single corner, which is consistent to the topological corner charge $Q^{(4)}_c=\frac{1}{4}$. It is same for the case of $C_2$ symmetric lattice (see Methods).

	In our experiment, we fabricate the lattice with parameters as $d_a^x=d_a^y$ varying from 13 to $13.4\ \mu$m with step of $0.1\ \mu$m, and $d_b^x=d_b^y$ fixed as $11\ \mu$m. We inject the single photon into lattices from one of the lattice corners, the single-photon is not localized in the excited corner as previous results shown in Fig.~\ref{f2} and the probability distribution occupies all four corners and the sites in the lattice boundaries. We confirm and demonstrate the result with four kind lattices with different parameter of $t_a/t_b$ varying from 0.65 to 0.70 and evolution distances as shown in Fig.~\ref{f4}(c). If we investigate the four-corner output probability distribution as the sub-space of the result, we can find that the finite gap effect can provide a way to obtain single-photon superposition state in a lattice.

	\paragraph*{Robustness.} The observed corner states are 2D analogy of Jackiw-Rebbi solitons and robust against disorders and perturbations as long as the rotation symmetry and time-reversal symmetry are preserved. From the perspective of topological band theory, for the $C_2$ lattice, the corner states are mid-gap states and separated from the bulk states. Therefore, as long as the disorders does not close the band gaps and breaking the symmetries, the corner states remain unchanged. For $C_4$ symmetric case, the corner states are not in the band gap, however, they are decoupled from the bulk bands and can be excited individually by photons injection on four rods simultaneously.

	Meanwhile, to demonstrate the performance of topological corner modes, we fabricate six samples and each sample also is realized in different evolution distance varying from 10 to 30 mm in step of 5 mm. It means that thirty lattices are fabricated and measured in our experiment. Though all the parameters of the system and the fabrication environment during femtosecond laser direct-write process have been locked and optimized, it is still inevitable that the disorder induced by fabrication exists in each lattice. Nevertheless, the performance of corner states in the experimental result almost have no change, implying the distinguishing robust feature of topological phase.

	Besides, we also introduce the disorder into the pure lattice and demonstrate the robustness of the topological phase. In our experiment, we set the $\eta=\varDelta d/\bar{d}$ as disorder level, where $\varDelta d$ is the introduced disorder separation distance and $\bar{d}$ is the averaged separation distance of the pure lattice. As shown in Fig.~\ref{s_robustness}, the photon is still well confined in the excited sites, implying the robustness of the explored second-order topological phase. We set $d_a^x = d_a^y = 14$ mm, $d_b^x  = d_b^yc  = 18$ mm for $C_4$ symmetric trivial lattice and $d_a^y = d_b^y = 16$ mm, $d_a^x = 14$ mm, $d_b^x  = 18$ mm for $C_2$ symmetric trivial lattice. The photons can not be confined in the corners and diffuse into the whole lattice, as shown in Fig.~\ref{s_robustness}.\\

	\subsection*{Conclusion and discussion}
	In summary, we present direct observation of higher-order topological phases and the topological protection to quantum superposition and entanglement in two-dimensional photonic lattice fabricated with femtosecond laser direct writing technique. We demonstrate the different properties of second-order corner modes in the photonic lattice with $C_4$ and $C_2$ symmetry.

	We employ single-photon superposition state to identify the corner states in the meanwhile show a topological protection mechanism isolating the quantum entanglement from diffusion-induced decoherence, which provides a promising and valuable resource for quantum information processing.

	Our implementation of SOTI at single-photon level can facilitate the way for studying lower dimensional topological localized states with controllable photon numbers in the quantum regime. Unlike recent works on $C_6$~\cite{HOTI10} and $C_3$~\cite{C3} symmetric waveguide arrays, we exhaust the predictions of photonic SOTI in the remaining $C_4$ and $C_2$ symmetric systems, and point out that the corner state can be embedded in the bulk states while being decoupled and can be excited individually from them. Our results extend the protection mechanism of topological phases into quantum regime, and demonstrate the compatibility of quantum photonics and topological phase. The demonstrated key elements, including integrated structures, higher-order topological crystalline insulator, built-in superposition state generation, and heralded single-photon imaging, can enrich the field of quantum topological photonics.\\

	\paragraph*{Acknowledgments}
	\noindent The authors thank Roberto Osellame and Jian-Wei Pan for helpful discussions. This research is supported by National Key R\&D Program of China (2017YFA0303700, 2017YFA0303702,2017YFA0304203, 2018YFA0306200, and 2019YFA0308700), National Natural Science Foundation of China (11690033, 61734005, 11761141014, 11604392, 11625418, 11890700, and 51732006), Science and Technology Commission of Shanghai Municipality (17JC1400403), Shanghai Municipal Education Commission (2017-01-07-00-02-E00049), IRT\_17R70, 1331KSC and 111 Project (D18001), China Postdoctoral Science Foundation Funded Project (2019M661784). X.-M.J. acknowledges additional support from a Shanghai talent program.\\

	\subsection*{Methods}

	\paragraph*{Topological classification and bulk polarization:}
	The time-reversal symmetry of waveguide lattice leads to a vanishing Berry curvature and a zero-Chern number~\cite{SSH2,SSH3}. However, the extra $C_2$ and $C_4$ rotation symmetries will put topologically non-trivial constrains on the eigenfunctions, forming the topological crystalline insulators (TCIs)~\cite{TCI1}. In our cases, the TCIs can be classified by the 2D bulk polarization~\cite{POL1,POL2} defined as follows,
	\begin{equation}
		P_i=-\frac{1}{(2\pi)^2}\int_{BZ} d^2\bm{k}\mathrm{Tr} [\hat{{\cal A}}_i] , \quad i=x,y
	\end{equation}
	where $BZ$ presents the first Brillouin zone, $(\hat{{\cal A}}_i)_{mn}(\textbf{k})=\mathrm{i}\bra{u_m(\textbf{k})}\partial_{k_i}\ket{u_n(\textbf{k})}$, with $m$, $n$ run over all bands below the gap, $\ket{u_m(\textbf{k})}$ is the periodic part of the eigenfunction for the $m$th band. The 2D polarization is simply related to the 2D Zak phase via $\theta_i=2\pi P_i$ for $i=x,y$.

	Besides, the value of the 2D bulk polarization is equal to position of the Wannier center~\cite{WC}. Due to the two mirror symmetries in the waveguide array, the bulk polarization is quantized and the Wannier center is restricted at the maximal Wycoff position of the unit cell. The Wannier center can be applied to investigate the topological classes. If the Wannier center is restricted at the center of the unit cell, namely $(P_x, P_y)=(0, 0)$, the system is adiabatically connected to the atomic insulators which is topologically trivial insulators. Nevertheless, if the Wannier center is restricted to the center of the edge of the unit cell, namely $(P_x, P_y)=(0, \frac{1}{2})$ or $(P_x, P_y)=(\frac{1}{2}, 0)$, it corresponds to the first-order topological insulator with 1D edge states. When $(P_x, P_y)=(\frac{1}{2}, \frac{1}{2})$ which means that the Wannier center is located at the corner of the unit cell, the system is a second-order topological insulator with both 1D edge states and 0D corner states. The topological edge states and corner states are the Jackiw-Rebbi solitons for the two topologically distinct bulks and edges respectively~\cite{soliton}.

	For $C_4$ symmetric array, we have $P_x=P_y$ and since $P_x$ can be either $\frac{1}{2}$ or $0$, the bulk polarization forms a $\mathbb{Z}_2$ topological index of the system. However, this is not case for $C_2$ symmetric array because $P_x$ and $P_y$ can be independent to each other. Therefore the bulk polarization forms a $\mathbb{Z}_2\times\mathbb{Z}_2$ topological index for $C_2$ symmetric system. For our cases, when $t^i_a<t^i_b$ for $i=x,y$, we have $(P_x, P_y)=(\frac{1}{2}, \frac{1}{2})$ which implies that the waveguide array is in the SOTI phase.

	\paragraph*{Determining the values of topological indices:} To further determine the values of topological indices, we need to investigate the coupling configurations of our lattice. If we consider two $n-$fold rotation symmetric topological crystalline insulator, the sum of them is also a TCI with the symmetry being the sum of previous $n-$ fold rotation symmetry. This property ensures a free Abelian additive structure of the classification of TCIs and therefore we can choose a set of primitive systems to generate all TCIs up to stable equivalence~\cite{TCI1}. We define these primitive systems as the primitive generators which satisfy certain rotation symmetry.

	For our $C_4$ symmetric lattice, the primitive generator is $h^{(4)}_{1b}$ according to the algebraic method~\cite{TCI1, TCI2}. The previous 2D bulk polarization can be directly obtained from the topological indices as
	\begin{equation}
		P^{(4)}_x=P^{(4)}_y=\frac{1}{2}[X_1],
	\end{equation}
	and it is defined modulo 1. The corner states arise due to the filling anomaly: the mismatch between the $C_4$-symmetry and conservation of the number of photonic eigenstates and we can define a topological corner index as follows
	\begin{equation}
		Q^{(4)}_c=\frac{1}{4}([X_1]+2[M^{(4)}_1]+3[M^{(4)}_2]),
	\end{equation}
	and it is also defined modulo 1. In $C_4$ symmetric lattice, for the non-trivial case, we have $[X_1]=-1$, $[M_1]=1$ and $[M_2]=0$. Therefore the bulk polarization is $P^{(4)}_x=P^{(4)}_y=\frac{1}{2}$ which is consist with previous calculations and the topological corner index is $Q^{(4)}_c=\frac{1}{4}$, indicating $\frac{1}{4}$ fractionalized corner states at each of four corners.

	For $C_2$ symmetric lattice, the primitive generator of our waveguide array is also $h^{(4)}_{1b}$. The corresponding 2D bulk polarization is now calculated by
	\begin{equation}
		P^{(2)}_x=-\frac{1}{2}([Y_1]+[M^{(2)}_1]),
		P^{(2)}_y=-\frac{1}{2}([X_1]+[M^{(2)}_1]),
	\end{equation}
	and the topological corner index is defined as
	\begin{equation}
		Q^{(2)}_c=\frac{1}{4}(-[X_1]-[Y_1]+[M^{(2)}_1]).
	\end{equation}
	In $C_2$ symmetric lattice, for the non-trivial case, we have $[X_1]=-1$, $[Y_1]=-1$ and $[M^{(2)}_1]=0$. Therefore the bulk polarization is $P^{(2)}_x=P^{(2)}_y=\frac{1}{2}$ which is consist with previous calculations. However, in this case, the topological corner index is $Q^{(2)}_c=\frac{1}{2}$, indicating $\frac{1}{2}$ fractionalized corner states at each of four corners which is different from the $C_4$ symmetric lattice.

	\paragraph*{Fabrication and measurement of the lattices on a photonic chip:} We fabricate the samples in alkaline earth boro-aluminosilicate glass substrate (refractive index $n_0=1.514$ for the writing laser at a wavelength of 513 nm) using the laser system operating at a repetition rate of 1 MHz and a pulse duration of 290 fs. The light is focused inside the sample with a 50X microscope objective (NA=0.50) after being reshaped with a spatial light modulator. We continuously move the substrates using a high-precision three-axis translation stage with a constant velocity of 10 mm/s to create the lattices by the laser-induced refractive index increase.

	In the experiment, we inject the photons into the input waveguides in the photonic chip using a 20X objective lens. After a total propagation distance through the lattice structures, the outgoing photons are first collimated with a 10X microscope objective, then detected and analyzed by a combination of wave plates and polarizers.

	\paragraph*{The generation and imaging of the heralded single-photon state:} The single-photon source with the wavelength of 810 nm are generated from periodically-poled KTP (PPKTP) crystal via type-II spontaneous parametric down conversion. The generated photon pairs are separated to two components, horizontal and vertical polarization, after a long-pass filter and a polarized beam splitter (PBS). One should notice that the measured patterns would come from the thermal-state light rather than single photons if we inject only one polarized photon into the lattices without external trigger. Therefore, we inject the horizontally polarized photon into the lattices, while the vertically polarized photon acts as the trigger for heralding the horizontally polarized photons out from the lattices with a time slot of 10 ns. The measured second-order anti-correlation parameter is 0.026$\pm$0.003, implying that a single photon is well preserved in the corner states. We capture each evolution result using the ICCD camera after accumulating in the external mode for 600s.

	\paragraph*{The 3D 1$\times$4 photonic coupler:}
	The Hamiltonian of the 3D 1$\times$4 photonic coupler could be written as
	\begin{equation}
		H=
		\left(
		\begin{array}{ccccc}
			0 & c & c & c & c\\
			c & 0 & 0 & 0 & 0\\
			c & 0 & 0 & 0 & 0\\
			c & 0 & 0 & 0 & 0\\
			c & 0 & 0 & 0 & 0\\
		\end{array}
		\right),
		\label{Hc}
	\end{equation}
	where the entry waveguide is labeled as 1 and the other four waveguides are labeled as 2 to 5 respectively, $c$ is the coupling strength. In our experiment, the distance of 1$\times$4 photonic coupler is set as $L = \frac{\pi}{4c}$. According to the evolution operator $U=e^{-iHL}$, we can obtain
	\begin{equation}
		U=
		-\frac{1}{4}\left(
		\begin{array}{ccccc}
			0 & 2i & 2i & 2i & 2i\\
			2i & -3 & 1 & 1 & 1\\
			2i & 1 & -3 & 1 & 1\\
			2i & 1 & 1 & -3 & 1\\
			2i & 1 & 1 & 1 & -3\\
		\end{array}
		\right).
		\label{coupler}
	\end{equation}
	When we inject the single photon in to the entry waveguide, then $\ket{\psi_{in}}=[1\ 0\ 0\ 0\ 0]^T$. According to $\ket{\psi_{out}}=U\ket{\psi_{in}}$, we obtain $\ket{\psi_{out}}=\frac{i}{2}[0\ 1\ 1\ 1\ 1]^T$, implying that the probability and phase of single photon in four waveguides are uniform.

	\paragraph*{The band gap of $C_2$ symmetric lattice:}
	For $|t^i_a/t^i_b| \ll 1$ in $C_2$ symmetric lattices, the band gap is large and the corner states are well separated from the bulk states. The real space field distributions of corner states decay away from the corner position in forms of $e^{-m\delta x}$ and $e^{-m\delta y}$ along $x$ and $y$ directions where $m$ stand for the mass terms in the effective Hamiltonian which is proportional to the size of the band gaps. $\delta x$ and $\delta y$ represent the distances from the corner positions along $x$ and $y$ directions respectively. Therefore, in this case, the photons are well confined at four corners. On contrary, when $|t^i_a/t^i_b|$ is gradually increased to 1, the size of band gaps reduced and the localization of corner states are weakened, leading to field distributions in the edge and bulk of the lattices.

	\clearpage


\section{Introduction}
The negatively charged nitrogen-vacancy (NV$^-$ or just NV) center in diamond has opened up a unique opportunity
for nanoscale nuclear magnetic resonance (NMR) spectroscopy and imaging.~\cite{MKS+13,SSP+13,MKC+14,HSR+15,DPL+15,LSU+16,APN+17}
NV's performance as a room temperature quantum sensor for nanoNMR is crucially dependent on
its depth from the diamond surface ($d_{\mathrm{NV}}$), as the NV center closer to the surface exhibits a stronger interaction with NMR analyte placed on the surface.
Furthermore, near-surface NV centers lead to efficient couplings with light,
which also adds impetus to other applications such as nanophotonics, plasmonics, and quantum communication.~\cite{TMO+13,SMZ+16,AEV+18,AHWZ18,WEH18,BGF+19}

There are several approaches to creating near-surface NV centers.
Most commonly used is low-energy N$^+$ ion implantation, typically with the acceleration energy less than 5~keV.~\cite{PNJ+10,OPC+12,dOAW+17,FBH+18,SDS+19}
Other methods such as nitrogen $\delta$-doping by chemical vapor deposition (CVD) and precise surface etching by plasma or thermal oxidation have also been demonstrated.~\cite{OHB+12,ORW+13,MDD+14,LPMD14,KMS+14,CGO+15,dOMW+15,ZZW+17}
Of particular interest to these platforms is the relation between $d_{\mathrm{NV}}$ and $T_2$, the coherence time of NV's electronic spin.
Longer $T_2$ is required to achieve better magnetic sensitivity, whereas near-surface NV centers tend to suffer short $T_2$.

In Ref.~\onlinecite{ISS+17}, some of the present authors reported that
near-surface NV centers can also be created by comparatively high energy (10~keV) N$^+$ ion implantation,
when combined with a SiO$_2$ screening mask deposited on the diamond surface.
The defining feature of this method is that the region of the highest N$^+$ density is located at the surface.
We showed this by Monte Carlo simulations, and experimentally confirmed the creation of NV centers
(see Sec.~\ref{sec_sample} for a more detailed description of our method and the result of Monte Carlo simulations).
The resulting N$^+$ distribution in our method is in contrast with the case for the standard (maskless) low-energy ion implantation,
in which the implanted N$^+$ ions are approximately Gaussian-distributed inside of the diamond.
An experimental confirmation of the distribution of the NV centers along the depth direction (depth profile),
as well as the relation between $T_2$ and $d_{\mathrm{NV}}$, are highly desired but have not yet been reported so far.
The present work addresses these issues.
Specifically, we focus on single NV centers and determine $d_{\mathrm{NV}}$ from NMR of proton ensembles on the surface, with direct application to nanoNMR in mind.~\cite{AS18,SIA18}

The rest of this paper is organized as follows.
In Sec.~\ref{sec_density}, we determine the density and yield of NV centers created by our method.
Section~\ref{sec_t2} presents the main result, the relation between $T_2$ and $d_{\mathrm{NV}}$,
and we analyze the depth profile of the NV centers by assembling the data from as many as 141 NV centers.
In Sec.~\ref{sec_noise}, we perform noise spectroscopy to investigate the source of the noise that limits $T_2$ in our sample.
The diamond surface is examined by atomic force microscopy (AFM) and X-ray photoelectron spectroscopy (XPS) in Sec.~\ref{sec_surface}.
We conclude in Sec.~\ref{sec_conclusion}.
To avoid digressing from the main subject, technical aspects and supporting data are presented in Appendices for interested readers.
They include the sample preparation procedure and the Monte Carlo simulations of N$^+$ ion implantation with screening masks (Sec.~\ref{sec_sample}),
our experimental setup and microwave pulse sequences (Sec.~\ref{sec_pulse}),
and the result from double electron--electron resonance (DEER) measurements (Sec.~\ref{sec_deer}).
Although many of them are now becoming common tools among the researchers working on NV-based quantum sensing, it will still be beneficial to provide concise summaries.
We also refer to Ref.~\onlinecite{MSI+20} as a resource for further details about our experimental setup and techniques extensively used in the present work.

\section{Density and yield\label{sec_density}}
Figures~\ref{fig1}(a--d) show fluorescence images (20~$\times$~20~$\mu$m$^2$) of diamond surfaces,
where the SiO$_2$ layers of thickness $t$ = 52.3, 57.6, 64.1, and 69.1~nm were deposited prior to N$^+$ ion implantation.
\begin{figure}
\begin{center}
\includegraphics{fig1.eps}
\caption{(a--d) Fluorescence images (20~$\times$~20~$\mu$m$^2$) of diamond surfaces under the laser power of 0.6~mW.
$t$: thickness of the SiO$_2$ layer.
kcps: kilo-counts per second.
(e) NV density as a function of $t$.
(f) Yield as a function of $t$.
\label{fig1}}
\end{center}
\end{figure}
It should be noted that these SiO$_2$ layers had been removed when these images were taken (see Sec.~\ref{sec_sample}).
Each image shows multiple bright spots, the number of which is decreasing as $t$ increases.
To confirm the presence of single NV centers, continuous wave optically detected magnetic resonance (CW ODMR) was performed on individual bright spots.
For the areas with $t$ = 64.1 and 69.1~nm, all the bright spots in the searched area were examined and the number of ODMR-active points were counted.
The areal density was then calculated by dividing it by the searched area.
For the $t$ = 52.3 and 57.6~nm areas, where many more bright spots exist, we set a threshold count rate, above which we judged as an emission from a single NV center,
and the number of above-threshold spots were counted by a computer program.
The threshold count rates were 2.0~$\times$~10$^4$~kcps for $t$ = 52.3~nm and 1.7~$\times$10$^4$~kcps for $t$ = 57.6~nm,
which were the minimum photon count rates among a subset of the bright spots on which CW ODMR was performed and the signals were observed.
The result is summarized in Fig.~\ref{fig1}(e).

We also calculate the yield (efficiency of conversion from N$^+$ ions into NV$^-$ centers) as shown in Fig.~\ref{fig1}(f).
It is given by the ratio of the areal density divided by the implanted N$^+$ density deduced from Monte Carlo simulations.
The low yield of 0.1--0.4\% is consistent with our previous observation,~\cite{ISS+17}
and similar to the values obtained in low energy ion implantation ($\approx$ 3~keV).~\cite{PNJ+10}
The cause of the latter case is attributed to the lack of vacancies to pair up with nitrogen atoms to form NV centers,
but it has been shown that the yield can be improved by carefully adjusting annealing, chemical or plasma etching conditions.~\cite{OSP+13,YUW+13,AHA+14}
This aspect has not been optimized yet in our case, and can be improved in the future.

\section{Coherence time and depth\label{sec_t2}}
Figure~\ref{fig2}(a) shows an example of Hahn echo decay [see Fig.~\ref{fig7}(a) for the pulse sequence].~\cite{H50}
\begin{figure}
\begin{center}
\includegraphics{fig2.eps}
\caption{(a) Hahn echo decay.
(b) NMR spectrum.
(c) Relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$.
The dashed line is the measurement limit as determined by Eq.~(\ref{eq_t2limit}), whereas the dash--dot line indicates $T_2$ = 3~$\mu$s.
They define the region (unpainted) that is ``accessible'' in our experiments.
See the main text for detail.
The filled circle indicates the NV center with which the Hahn echo decay in (a), the NMR spectrum in (b), and the noise spectroscopy data in Fig.~\ref{fig4} were taken.
The filled triangle indicates the NV center with which the DEER data in Fig.~\ref{fig8} were taken.
\label{fig2}}
\end{center}
\end{figure}
The static magnetic field $B_0$ = 23.2~mT was applied parallel to the NV symmetry axis,
and the microwave frequency was set at $\omega_{0,-1}/2\pi$ = 2218.2~MHz, corresponding to the $m_S = 0 \leftrightarrow -1$ transition.
$P_0$ is the probability of the final spin state being $m_S$ = 0, which we calculate from calibrated photon counts.~\cite{MSI+20}
We describe a general coherence decay curve as
\begin{equation}
P_0 = \frac{1}{2} + \frac{ \mathcal{C} }{2} \exp \left[ -\left( \frac{ t_{\mathrm{tot}} }{T_2} \right)^p \right].
\label{eq_t2}
\end{equation}
The form of $\mathcal{C}$ depends on a type of measurements, $t_{\mathrm{tot}}$ is the total time of the pulse sequence applied.
Fitting parameters are $T_2$ and $p$ (stretched exponent).
In the case of Hahn echo, $\mathcal{C}_{\mathrm{echo}}$ = 1 and $t_{\mathrm{tot}}$ = 2$\tau$.
The solid line in Fig.~\ref{fig2}(a) is the fit, giving $T_{2,\mathrm{echo}}$ = 27.1~$\mu$s and $p$ = 0.98.

We then apply a multipulse sequence to determine $d_{\mathrm{NV}}$ from NMR of protons on the surface.
We use XY4, XY8, and XY16 sequences (see Sec.~\ref{sec_pulse} for the definitions).~\cite{GBC90}
When the total number of $\pi$ pulses in a sequence is $N$, it is denoted as XY$k$-$N$ ($k$ = 4, 8, or 16 and $N$ is an integer multiple of $k$).
The decay is given by setting $t_{\mathrm{tot}}$ = $N\tau$, and the NMR spectrum is described by~\cite{PDC+16}
\begin{eqnarray}
& & \mathcal{C}_{\mathrm{NMR}} = \\
& & \exp \left\{ - 2 \left( \frac{ \gamma_{\mathrm{e}} B_{\mathrm{rms}} N\tau }{ \pi } \right)^2
\mathrm{sinc}^2 \left[ \frac{N \tau}{2} \left( \omega_{\mathrm{n}} - \frac{\pi}{\tau} \right) \right] \right\}. \nonumber
\end{eqnarray}
Here, $\gamma_{\mathrm{e}}/2\pi$ = 28~MHz/mT is the gyromagnetic ratio of the NV spin, and
$\omega_{\mathrm{n}}/2\pi = \gamma_{\mathrm{h}} B_0/2\pi$ is the proton NMR frequency
with $\gamma_{\mathrm{h}}/2\pi$ = 42.577~kHz/mT the gyromagnetic ratio of the $^1$H nuclear spin.
$B_{\mathrm{rms}}$ is the root-mean-square nuclear dipolar magnetic field integrated over the semi-infinite volume (i.e., outside of the diamond),
and can be calculated analytically as~\cite{SSP+13,LPMD14,RDO+14,PDC+16}
\begin{equation}
B_{\mathrm{rms}} = \frac{\mu_0 \hbar \gamma_{\mathrm{h}} }{4 \pi} \sqrt{ \frac{ 5 \pi \rho}{ 96 \, d_{\mathrm{NV}}^3 } }
\end{equation}
with $\mu_0$ the permeability of vacuum, $\hbar$ the reduced Planck constant, and
$\rho$ the proton density of the immersion oil on the surface, known to be 6~$\times$~10$^{28}$~m$^{-3}$ in the present case (Olympus IMMOIL-F30CC).

Figure~\ref{fig2}(b) is the NMR spectrum, taken with the NV center of Fig.~\ref{fig2}(a) and the XY16-128 sequence.
The stretched exponential decay component is subtracted (normalized $P_0$, or $P_{0,\mathrm{norm}}$),
and the horizontal axis is $(2\tau)^{-1}$ in the frequency unit.
The NMR dip appears at 0.988~MHz, agreeing with the calculated $\omega_{\mathrm{n}}/2\pi$.
From the fit, $d_{\mathrm{NV}}$ is determined to be 17.4~nm.
(Consult Fig.~10 of Ref.~\onlinecite{MSI+20} for another example of the NMR spectrum measured with a different, shallower NV center in the same sample.)

We repeat the measurements on different NV centers, and the relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$,
obtained from 43 NV centers in total, is plotted in Fig.~\ref{fig2}(c).
[The point with the filled circle corresponds to the NV center examined in Figs.~\ref{fig2}(a,b).]
We note again that, although the NV centers are grouped by $t$ (SiO$_2$ thickness), the SiO$_2$ layer had been removed when these measurements were performed.
A general trend is clearly discerned: the deeper the NV center, the longer the coherence,
as observed in previous works using different methods to create NV centers.~\cite{dOAW+17,FBH+18,SDS+19}
For instance, we find the observed $T_{2,\mathrm{echo}}$--$d_{\mathrm{NV}}$ relation is quantitatively close to that obtained in Ref.~\onlinecite{SDS+19} (without high temperature annealing).
It must be emphasized, however, that in our case a large fraction of NV centers exhibited short $T_2$ insufficient to observe clear NMR spectra.
Consequently, we could not determine their depths, and such NV centers are not shown in Fig.~\ref{fig2}(c).
We will include them when estimating the depth profile,
but before doing so we discuss an experimental condition to detect NMR and set the region in the $T_{2,\mathrm{echo}}$--$d_{\mathrm{NV}}$ space
where the experimental observation can be made.

The condition is simply that the signal $\mathcal{S}$ be larger than the noise $\mathcal{N}$.
The NMR signal arises from the component of $P_0$ after subtracting the coherence decay,
therefore $\mathcal{S} = \frac{1}{2} (1-\mathcal{C}_{\mathrm{NMR}}) e^{-(N\tau/T_2)^p}$.
We set $p$ = 1 for simplicity (also consistent with experimental observation) and $\tau = \pi/\omega_{\mathrm{n}}$ = 506.2~ns, corresponding to the resonance condition.
Experimentally, the number of $\pi$ pulses $N$ is determined from the balance between the signal strength and the coherence time,
and we usually found $N\tau/T_{2,\mathrm{echo}} \approx 2$ to be an appropriate operating condition.
For instance, the condition for Fig.~\ref{fig2}(b) is $N\tau/T_{2,\mathrm{echo}}$ = 128~$\times$~506.2~ns/27.1~$\mu$s $\approx$ 2.4.
Here, $T_{2,\mathrm{echo}}$ is used, as we do not know $T_2$ under dynamical decoupling in advance.

The measurement noise $\mathcal{N}$ arises from the Poisson distribution of the photon counts $n$, the variance of which is given by $\sqrt{n}$.
The signal photon count $n_{\mathrm{s}}$ is defined as $n_{\mathrm{s}} = n_1 + (n_0 - n_1) P_0$,
where $n_0$ and $n_1$ are the photon counts from the $m_S$ = 0 and $-$1 states, respectively, and are mutually related by $n_1 = (1-c) n_0$, with $c$ the contrast.
We obtain $P_0 = (n_{\mathrm{s}} - n_1)/(c n_0) $
and the fluctuation of $n_{\mathrm{s}}$, given by $\sqrt{n_{\mathrm{s}}}$, leads to $\mathcal{N} = \sqrt{ n_{\mathrm{s}} }/(c n_0) \approx (c \sqrt{n_0})^{-1}$.
Combining these, we write the condition $\mathcal{S}/\mathcal{N} = 1$ as
\begin{equation}
T_{2,\mathrm{echo}} = \frac{ \pi }{ 2\sqrt{2} \gamma_e B_{\mathrm{rms}} } \sqrt{-\ln \left( 1- \frac{2 e^2}{ c \sqrt{n_0} } \right) }.
\label{eq_t2limit}
\end{equation}
Typical values for $c$ and $n_0$ in the present experimental condition are $c \approx$ 0.2 and $n_0 \approx 5 \times 10^5$ counts.
Substituting these values, we draw the dashed line in Fig.~\ref{fig2}(c), below which the depth cannot be determined.
In addition to this constraint, short $T_{2,\mathrm{echo}}$ limits the number of pulses we can apply.
For instance, to satisfy $N\tau/T_{2,\mathrm{echo}} \approx 2$ with $N$ = 8, $T_{2,\mathrm{echo}}$ needs to be longer than 2~$\mu$s.
In fact, from NV centers with $T_{2,\mathrm{echo}}$ shorter than 3~$\mu$s we were not able to observe clear NMR signals
(and the extension of $T_2$ by dynamical decoupling is often not effective),
even though they are likely to be located close to the surface.
This motivates us to draw another bound $T_{2,\mathrm{echo}}$ = 3~$\mu$s in Fig.~\ref{fig2}(c).
Although the argument here is manifestly crude, the defined ``accessible'' region (unpainted) reasonably coincides with the experimental observation.

From the dataset in Fig.~\ref{fig2}(c), it is straightforward to produce histograms
showing the number of NV centers (occurrence) as a function of $d_{\mathrm{NV}}$ [Figs.~\ref{fig3}(a--d)].
\begin{figure}
\begin{center}
\includegraphics{fig3.eps}
\caption{(a--d) Histograms showing the distributions of the 43 {\it depth-determined} NV centers.
(e--h) NV densities estimated from the experiments ($\bigcirc$), 98 {\it depth-undetermined} NV centers included, and the simulations (solid lines).
See the main text for detail.
\label{fig3}}
\end{center}
\end{figure}
Ideally, these histograms are to be directly compared with the depth profile from Monte Carlo simulations (Fig.~\ref{fig6}).
However, as mentioned above, the fact that the data points in Fig.~\ref{fig2}(c) are limited to a certain region suggests
that they alone are not sufficient to deduce a firm conclusion about the depth profile.
There are indeed another 98 NV centers we measured but $d_{\mathrm{NV}}$ were not determined.
This number breaks down into 44 ($t$ = 52.3~nm), 33 (57.6~nm), 10 (64.1~nm), and 11 (69.1~nm),
not negligible compared with the numbers of the depth-determined NV centers, 13 (52.3~nm), 12 (57.6~nm), 9 (64.1~nm), and 9 (69.1~nm).

To include the former group in the histograms, we assume, judging from the trend observed in Fig.~\ref{fig2}(c), that they are located at $d_{\mathrm{NV}} \leq$ 5~nm.
In Figs.~\ref{fig3}(e--h), the solid lines are the Monte Carlo simulations in Fig.~\ref{fig6} multiplied by the respective yields in Fig.~\ref{fig1}(f).
The circle points ($\bigcirc$) for $d_{\mathrm{NV}} >$ 5~nm are the histograms multiplied by the respective NV densities.
Similarly, those for $d_{\mathrm{NV}} \leq$ 5~nm are based on the numbers of the NV centers whose depths left undetermined, but divided by five to reflect the fact that
we are agnostic on which ``bin'' (the width of 1~nm) they are in.
Thus, the horizontal error bar for them is set as $d_{\mathrm{NV}}$ = 2.5~$\pm$ 2.5~nm (the other error bars are set as the bin width).
From Figs.~\ref{fig3}(e--h), it is interpreted, though not conclusive, that the profile is concentrated toward the surface,
unlike the standard ion implantation without a screening mask.

\section{Noise spectroscopy\label{sec_noise}}
In our sample, it is observed that the NV centers with short $T_{2,\mathrm{echo}}$ also tend to be short-lived, degrading $T_{2,\mathrm{echo}}$ during the experiments.
In one case, initially $T_{2,\mathrm{echo}}$ was 6~$\mu$s, but it reduced continuously to 2~$\mu$s in 5 hours, then fluctuates between 1 and 2~$\mu$s.
Another case shows that initial $T_{2,\mathrm{echo}}$ of 18~$\mu$s slowly degrades to 12~$\mu$s in the time scale of 80~hours.
These are likely to arise from the coupling with the surface, and we next examine the origin of the noise that affects the NV spin.
For this purpose, we would ideally hope to use shallower NV centers, but as mentioned above, they suffer degradation of $T_{2,\mathrm{echo}}$,
making a clean set of experiments difficult.
Nonetheless, the trend in Fig.~\ref{fig2}(c) indicates that the plotted NV centers share the same noise source.
We therefore chose to use the NV center with longest $T_{2,\mathrm{echo}}$ found in the $t$ = 52.3~nm area [the filled circle in Fig.~\ref{fig2}(c)].
For the noise spectroscopy, $B_0$ was fixed at 23.2~mT, where $T_2$ and $d_{\mathrm{NV}}$ had also been measured.

We first perform dynamical decoupling by increasing $N$ up to 512.
The decay curves are shown in Fig.~\ref{fig4}(a).
\begin{figure}
\begin{center}
\includegraphics{fig4.eps}
\caption{(a) Coherence decay curves under dynamical decoupling.
The $N$ = 1 data were taken by the Hahn echo (the time axis is given by 2$\tau$),
the $N$ = 8 data by XY8, and the $N \geq$ 16 data by XY16-$N$ (the time axis is given by $N\tau$).
The inset: $T_2$ as a function of $N$.
The horizontal bars indicate the 95\% confidence limits.
(b) Relaxation signals for the SQ ($\triangle$) and DQ ($\square$) channels.
(c) Noise spectrum $S_{\mathrm{e}}(\omega)$ deduced from the spectral decomposition of the dynamical decoupling data ($\bigcirc$),
the SQ relaxation rate $\Omega$ ($\triangle$), and the DQ relaxation rate $\gamma$ ($\square$).
All the measurements were performed at $B_0$ = 23.2~mT.
\label{fig4}}
\end{center}
\end{figure}
Although dynamical decoupling improves the coherence to some extent, the effect is not significant.
$T_2^{(N)}$, $T_2$ by $N$-pulse dynamical decoupling, shows a scaling and saturating behavior described by
\begin{equation}
T_2^{(N)} = T_2^{(1)} [ N_{\mathrm{sat}}^s + (N^s - N_{\mathrm{sat}}^s) e^{-N/N_{\mathrm{sat}}} ],
\end{equation}
and the fit gives $N_{\mathrm{sat}}$ = 105.2, $s$ = 0.10, and $T_2^{(\infty)}$ = $T_2^{(1)} N_{\mathrm{sat}}^{s}$ = 43.0~$\mu$s [the inset of Fig.~\ref{fig4}(a)].
$T_2^{(1)}$ corresponds to $T_{2,\mathrm{echo}}$, but is treated as a parameter here.
The scaling exponent $s$ of 2/3 is expected for an environment that fluctuates slowly compared with $T_2$.~\cite{BPB+12,RMU+15}
The small $s \ll$ 2/3 as well as the saturated $T_2$ observed here indicate the presence of broadband noise.

This observation is reinforced by spectral decomposition.
Under a few simplifying assumptions, the spectral component at $\omega = \pi/\tau$ is given by~\cite{BPB+12,RMU+15}
\begin{equation}
S_{\mathrm{DD}}(\omega) = -\frac{ \pi \ln(2P_0 -1)}{ N \tau }.
\end{equation}
One can compute the noise spectrum directly from the dynamical decoupling data.
Figure~\ref{fig4}(c) shows the spectral function $S_{\mathrm{e}}(\omega)$,
related to $S_{\mathrm{DD}}(\omega)$ by $S_{\mathrm{e}}(\omega) = 2S_{\mathrm{DD}}(\omega)/d_{\parallel}^2$,
where $d_{\parallel}$ = 0.35~Hz~cm/V is the electric dipole moment parallel to the NV symmetry axis.~\cite{MAJ17}
The reason for this conversion is to enable a comparison between the (effective) electric and magnetic field noises.
Note that in Fig~\ref{fig4}(c) the data points close to $(2P_0 -1)$ = 0 or 1 that are insensitive to the change in $\tau$ have been removed.
Further analysis usually invokes a Lorentzian fit~\cite{BPB+12,RMU+15}
\begin{equation}
S_{\mathrm{DD}}(\omega) = \frac{ \langle E_{\perp}^2 \rangle \tau_{\mathrm{c}} }{ \pi [1+(\omega \tau_{\mathrm{c}})^2] }
\end{equation}
or a sum of multiple Lorentzians,
where $\langle E_{\perp}^2 \rangle$ is the average coupling strength of the environment to the NV center,
and $\tau_{\mathrm{c}}$ is the correlation time of the environment.
In the present case, the spectrum is almost flat (typically 1--2$\times$10$^{10}$~V$^2$/m$^2$Hz)
up to $\omega/2\pi \approx$ 10$^7$~Hz in the double logarithmic plot,
suggesting that the noise is white in the probed frequency range ($\omega \tau_{ \mathrm{c}} \ll 1$).
This makes fitting to a Lorentzian curve arbitrary, and we can only infer $\tau_{ \mathrm{c}}$ to be shorter than 100~ns.
Still, the short $\tau_{ \mathrm{c}} \ll T_2$ is consistent with the small $s$.

We also probe the noise by relaxation measurements through single quantum (SQ) and double quantum (DQ) channels
as shown in Fig.~\ref{fig4}(b) [see Fig.~\ref{fig7}(c) for the pulse sequence].~\cite{MAJ17}
The relaxation times of the two channels are almost the same,
$T_{1,\mathrm{SQ}}$ = 1.06~ms and $T_{1,\mathrm{DQ}}$ = 1.16~ms,
consistent with the previous observations at similar magnetic fields.~\cite{MAJ17,SDS+19,GCK20}
We can then deduce the SQ and DQ relaxation rates, $\Omega$ and $\gamma$,
using the relations $T_{1,\mathrm{SQ}} = (3 \Omega)^{-1}$ and $T_{1,\mathrm{DQ}} = (\Omega+2\gamma)^{-1}$.~\cite{MAJ17}
We obtain $\Omega$ = 316~Hz, corresponding to the magnetic noise probed at $\omega_{0,-1}/2\pi$ = 2218.2~MHz,
and $\gamma$ = 272~Hz, corresponding to the electric noise probed at $\omega_{-1,1}/2\pi$ = 1301.5~MHz (the separation of the $m_S$ = $-$1 and 1 states at 23.2~mT).
$\Omega$ and $\gamma$ are converted into the spectral components as
$2 \Omega/d_{\parallel}^2$ = 5.16$\times$10$^7$~V$^2$/m$^2$Hz and
$\gamma/d_{\perp}^2$ = 9.41$\times$10$^{3}$~V$^2$/m$^2$Hz, respectively,
with $d_{\perp}$ = 17~Hz~cm/V the electric dipole moment perpendicular to the NV symmetry axis.~\cite{MAJ17,SDS+19}
The combined plot in Fig.~\ref{fig4}(c) certifies that in the wide frequency region the dominant contribution comes from the magnetic noise,
which is presumably due to surface defects, as has been argued in the case of low energy ion implantations.~\cite{RMU+15,MAJ17,SDS+19}

\section{Surface condition\label{sec_surface}}
We have found that our method of creating near-surface NV centers provides the $T_2$--$d_{\mathrm{NV}}$ relation comparable to other existing methods.
Yet, the SiO$_2$ deposition, in which an amorphous SiO$_2$ source (melting point $\approx$ 1700$^{\circ}$C) is heated by an electron beam
and sublimated SiO$_2$ is directed to the diamond substrate in vacuum,
is the process unique to ours, and it is of interest to examine an effect it might have on the diamond surface.
Here, we perform AFM and XPS on two areas, one in which SiO$_2$ was not deposited and ion implantation was blocked (denoted as $t$ = 0),
and the other in which ion implantation was performed through SiO$_2$ layers (denoted as $t \neq$ 0).
These measurements were performed after all the ODMR-based experiments were done (including DEER presented in Sec.~\ref{sec_deer}),
and therefore the $t$ = 0 area was subject to the same chemical and annealing processes (described in Sec.~\ref{sec_sample}) as the $t \neq$ 0 area.

Figure~\ref{fig5}(a) compares AFM topographic images (500~$\times$~500~nm$^2$) taken at the $t$ = 0 and $t \neq$ 0 areas.
\begin{figure}
\begin{center}
\includegraphics{fig5.eps}
\caption{(a) AFM topographic images (500~$\times$~500~nm$^2$) of the diamond surfaces.
(b) The normalized XPS spectra showing the carbon 1$s$ (left) and oxygen 1$s$ (right) peaks.
The spectra for $t \neq$ 0 are shifted by 0.1 for clarity.
\label{fig5}}
\end{center}
\end{figure}
To be more precise, the latter belongs to the $t$ = 52.3~nm area, but is not correlated with the area on which the ODMR experiments were carried out.
In these particular areas shown in Fig.~\ref{fig5}(a), the root-mean-square roughness, $R_{\mathrm{rms}}$, is 0.11~nm for $t$ = 0 and 0.14~nm for $t \neq$ 0.
The measurements were repeated on different locations, and $R_{\mathrm{rms}}$ averaged over the measured areas were 0.17~nm for both $t$ = 0 and $t \neq$ 0,
indicating that the SiO$_2$ deposition and subsequent ion implantation have not altered the surface condition significantly.
Nonetheless, the change in the surface morphology, with more `wrinkles' in the latter, is discernible.
It is not clear yet whether these wrinkles affect the coherence properties of the NV spins.

Figure~\ref{fig5}(b) shows XPS spectra of carbon 1$s$ and oxygen 1$s$.
The C 1$s$ peak at 285~eV arises primarily from $sp^3$ (in bulk),
and the asymmetric shape indicates the presence of other bonds such as
C--O, C=O ($>$ 285~eV), and $sp^2$ ($<$ 285~eV).~\cite{BDM+14,SDS+19}
The O 1$s$ peak at 532~eV has a broad linewidth, presumably composed of
C--O--C, C--C--H, and C=O ($<$ 532~eV).~\cite{BDM+14,SDS+19}
Despite the limited spectral resolution, we can conclude that the chemical composition of the surface is essentially the same between the $t$ = 0 and $t \neq$ 0 areas.

\section{Conclusion and outlook\label{sec_conclusion}}
To conclude, we created single nitrogen-vacancy (NV) centers by 10-keV N$^+$ ion implantation into diamond via SiO$_2$ layers of $t$ = 52.3, 57.6, 64.1, and 69.1~nm,
and examined the relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$.
The depth profiles incorporating the NV centers with $T_2 <$ 3~$\mu$s were found to be consistent with Monte Carlo simulations.
The effect of the surface on $T_2$ was examined by a combination of noise spectroscopy, surface topography, and XPS.

Although our method is found to give the $T_2$--$d_{\mathrm{NV}}$ relation comparable to other existing methods,
we also conceive that the surface undergoes a slight change in morphology due to SiO$_2$ deposition and $T_2$ becomes unstable and/or is shortened over time.
Recently, these aspects have been investigated in great detail by Sangtawesin {\it et al.},~\cite{SDS+19}
who emphasized the importance of keeping the surface purity and the quality of surface morphology in each step of sample preparation.
Longer coherence was achieved by combining annealing at 800$^{\circ}$C, at 1200$^{\circ}$C, and in oxygen atmosphere,
which are thought to remove residual vacancy complexes due to ion implantation and lead to highly ordered oxygen terminated surface.
While we did not adopt annealing at 1200$^{\circ}$C, this additional process is reported to improve $T_2$ several times longer (e.g., a few $\mu$s to ten $\mu$s)
and achieve $R_{\mathrm{rms}}$ of 0.06~nm (compared with our $R_{\mathrm{rms}}$ = 0.17~nm).
Given that without 1200$^{\circ}$C the $T_2$--$d_{\mathrm{NV}}$ relation in the present work and that in Ref.~\onlinecite{SDS+19} are similar,
we expect the process to provide an important clue for stabilizing and extending $T_2$ of NV center within 5~nm from the surface.

\section*{Acknowledgement}
K.S. was supported by the JSPS Grant-in-Aid for Research Fellowship for Young Scientists (DC1), Grant No.~JP17J05890.
T.T. was supported by the JSPS Grant-in-Aid for Scientific Research (KAKENHI) (B) Grant No.~20H02187 and (B) Grant No.~19H02617,
(S) Grant No.~16H06326
the JST CREST (JPMJCR1773), and MEXT Q‐LEAP (JPMXS0118068379).
K.M.I. was supported by the JSPS KAKENHI (S) Grant No.~26220602 and (B) Grant No.~19H02547,
the JST Development of Systems and Technologies for Advanced Measurement and Analysis (SENTAN),
and the Spintronics Research Network of Japan (Spin-RNJ).
The authors thank the technical staff of Central Service Facilities for Research at Keio University
for their assistance in XPS measurements.
The data that support the findings of the present work are available from the corresponding author upon reasonable request.