# Few-cycle excitation of atomic coherence: A closed-form analytical solution
beyond the rotating-wave approximation

Nazar Pyvovar  Bing Zeng  Lingze Duan<EMAIL_ADDRESS>Department of Physics
and Astronomy, The University of Alabama in Huntsville, Huntsville, AL 35899,
USA.

###### Abstract

Developing an analytical theory for atomic coherence driven by ultrashort
laster pulses has proved to be challenging due to the breakdown of the
rotating wave approximation (RWA). In this paper, we present an approximate,
closed-form solution to the Schrödinger equation that describes a two-level
atom under the excitation of a far-off-resonance, few-cycle pulse of arbitrary
shape without invoking the RWA. As an example of its applicability, an
analytical solution for Gaussian pulses is explicitly given. Comparisons with
numerical solutions validate the accuracy our solution within the scope of the
approximation. Finally, we outline an alternative approach that can lead to a
more accurate solution by capturing the nonlinear behaviors of the system.

††preprint: APS/123-QED

##  I Introduction

Quantum coherent control (QCC) is of great importance in fundamental physics
as well as a breadth of emerging applications [1]. With the advent of
femtosecond and attosecond light sources, control of atomic coherence using
ultrafast laser pulses with few optical cycles has attracted growing interest
in recent years [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Apart from
its significance in quantum theories, ultrafast QCC has profound implications
in practical applications. For example, in certain QCC schemes, using
ultrashort, broadband pulses allows the first electronic states of molecules
to be accessible and, at the same time, enables fast population transfer that
occurs well within the typical collision times [16, 17]. Few-cycle pulses can
also excite coherence on high-frequency transitions that enables efficient
generation of extreme ultraviolet (XUV) radiations [18, 12].

Studying ultrafast QCC in the few-cycle regime faces unique challenges. The
ultrashort pulse width invalidates the slowly-varying envelope approximation
(SVEA) [19] while the high peak field causes breakdown of the rotation-wave
approximation (RWA) [13]. As a result, the well-established theoretical
framework based on the optical Bloch equations and the area theorem ceases to
apply [20, 21]. Theoretical analysis has to rely on the Bloch equations or the
Schrödinger equation in their original forms without simplifications, which
are often highly nonlinear. This significantly increases the difficulty of
developing analytical theories. In most cases, numerical simulations have to
be used when dealing with few-cycle excitations [22, 23, 21, 24].

Meanwhile, there has been a continued effort to develop analytical theories
for atomic coherence driven by few-cycle pulses [2, 7, 8, 9, 10, 11, 12, 13,
14, 15]. A notably successful theory, proposed first by Rostovtsev et al.,
considers the coherence of a two-level atom under the excitation of a far-off-
resonance strong ultrashort pulse [9]. Through a perturbative scheme, the
model gives rise to a general solution of the Schrödinger equation without
invoking the RWA. The original solution, however, is not in closed form, and
the analysis of its features useful for practical applications still has to
rely largely on numerical computations. Several attempts have been made to
derive more explicit solutions under specific conditions [11, 25, 15]. In
particular, we have recently shown that a simple, closed-form analytical
solution of the Schrödinger equation can be obtained for few-cycle square
pulses [15].

In this paper, we expand our model to include arbitrary pulse shapes while
retaining an explicit, closed-form solution. Excitation by a few-cycle
Gaussian pulse is analyzed as an example of the general solution. The accuracy
of our solution is verified by comparing it to the exact numerical solution of
the general Schrödinger equation. Finally, an alternative approach is
suggested to simplify the theory, which can potentially lead to a more
accurate solution with a closer representation of the nonlinear behaviors of
the system and a broader scope of applicability.

##  II Closed-Form Solution

###  II.1 General Model

Our general model follows the theoretical framework described in Ref. [9]. A
quick outline is given below. We consider a two-level system (TLS) under the
influence of an electromagnetic field. The Hamiltonian of the system is

$\hat{H}=\hbar\omega_{c}\ket{c}\bra{c}-\mu\mathcal{E}(t)\ket{c}\bra{d}-[\mu\mathcal{E}(t)]^{*}\ket{d}\bra{c},$
(1)

where $\ket{c},\ket{d}$ are upper and lower levels, respectively, $\omega_{c}$
is the transition frequency, $\mathcal{E}(t)$ is the electric field, and $\mu$
is the dipole moment of the system. We are interested in the electric field of
the form $\mathcal{E}(t)=\Omega(t)\cos(\omega t+\phi)$, where $\Omega(t)$ is
the pulse envelope function and $\phi$ is carrier-envelope phase (CEP).

With this Hamiltonian, the equations of motion for the system are given by

$\displaystyle\dot{C}(t)=-i\Omega(t)\cos(\omega t+\phi)e^{i\omega_{c}t}D(t),$
(2a) $\displaystyle\dot{D}(t)=-i\Omega^{*}(t)\cos(\omega
t+\phi)e^{-i\omega_{c}t}C(t),$ (2b)

where $C(t)$ and $D(t)$ are the probabilities of $\ket{c}$ and $\ket{d}$,
respectively, i.e., $\ket{\Psi}=C(t)e^{-i\omega_{c}t}\ket{c}+D(t)\ket{d}$.

It proves useful to introduce the following quantity to simplify our equations
at this point

$\theta(t)=\int_{-\infty}^{t}\Omega(t^{\prime})\cos(\omega
t^{\prime}+\phi)e^{i\omega_{c}t^{\prime}}dt^{\prime}.$ (3)

With this definition, the equations (2a), (2b) become

$\displaystyle\dot{C}(t)=-i\dot{\theta}(t)D(t),$ (4a)
$\displaystyle\dot{D}(t)=-i\dot{\theta}^{*}(t)C(t).$ (4b)

In this paper, we are only interested in the case when electromagnetic field
exciting the system is an ultrashort pulse that is non-zero within a finite
time interval. For this reason, all quantities hereafter are considered only
in a finite time interval $t\in[0,\tau]$. In particular, $\Omega(t)$ is only
non-zero inside of this interval, and hence all the dynamics of the system
occur exclusively within the interval as well. With this in mind, we can
replace the lower integration limit in (3) with $0$.

By introducing the quantity $f=C/D$, the Schrödinger equation (4) can be
simplified to

$\dot{f}(t)=i\dot{\theta}^{*}(t)f^{2}(t)-i\dot{\theta}.$ (5)

The zeroth-order solution of (5) is obtained by neglecting the $f^{2}(t)$ term
[9], which yields

$f_{0}(t)=-i\theta(t).$ (6)

Now, define the first-order approximation $f_{1}(t)$ as the solution of (5)
that satisfies $(f_{1}(t)-f_{0}(t))^{2}\ll 1$. It is then straightforward to
show that $f_{1}(t)$ satisfies the equation

$\dot{f}_{1}(t)=2\theta(t)\dot{\theta}^{*}(t)f_{1}(t)+i\theta^{2}(t)\dot{\theta}^{*}(t)-i\dot{\theta}(t).$
(7)

In this way, we can define a sequence of approximations to (5). Specifically,
if the $k$th-order approximation $f_{k}$ is found, then the $(k+1)$th-order
approximation $f_{k+1}$ is the solution of $f$ that meets the condition
$(f_{k+1}-f_{k})^{2}\ll 1$, which satisfies the equation

$\dot{f}_{k+1}=2if_{k}(t)\dot{\theta}^{*}(t)f_{k+1}(t)-if_{k}^{2}(t)\dot{\theta}^{*}(t)-i\dot{\theta}(t).$
(8)

For far off-resonance, however, it has been shown that the first approximation
(7) is already quite accurate. Therefore, we focus on (7) in the current work.
The general solution to (7) has been found to be [15]

$f_{1}(t)=-\frac{i}{2}[\theta(t)+\int_{0}^{t}\dot{\theta}(t^{\prime})e^{\alpha(t^{\prime},t)}dt^{\prime}],$
(9)

where
$\alpha(t^{\prime},t)=2\int_{t^{\prime}}^{t}\theta(t^{\prime\prime})\dot{\theta}^{*}(t^{\prime\prime})dt^{\prime\prime}$.

###  II.2 Approximate Solution

Although the solution (9) is generic, the integral on the right-hand side
cannot be evaluated in closed-form for any realistic pulse shape.

To find a sufficiently accurate closed-form solution of (7), we can expand
$\alpha(t^{\prime},t)$ in powers of $\delta=t-t^{\prime}$ and approximate

$\alpha(t^{\prime},t)=2\int_{t^{\prime}}^{\delta+t^{\prime}}\theta\dot{\theta}^{*}dt^{\prime\prime}=2\theta(t)\dot{\theta}^{*}(t)\cdot[t-t^{\prime}],$
(10)

to the leading term in $\delta$.

Using the fact that all quantities considered here are assumed to be limited
to a finite time interval $t\in[0,\tau]$, we can further approximate the
leading term above as

$\alpha(t^{\prime},t)=2(t-t^{\prime})\frac{1}{\tau}\int_{0}^{\tau}\theta(t^{\prime\prime})\dot{\theta}^{*}(t^{\prime\prime})dt^{\prime\prime}=-i\alpha_{0}(t-t^{\prime}),$
(11)

where

$\alpha_{0}=\frac{2i}{\tau}\int_{0}^{\tau}\theta(t^{\prime\prime})\dot{\theta}^{*}(t^{\prime\prime})dt^{\prime\prime}$
(12)

is a constant (the average of $\alpha(t^{\prime},t)$), and $\tau$ is once
again the width of the few-cycle driving pulse. This allows us to simplify (9)
by noting

$\displaystyle\int_{-\infty}^{t}\dot{\theta}(t^{\prime})e^{\alpha(t^{\prime},t)}dt^{\prime}$
$\displaystyle=\int_{-\infty}^{t}\Omega(t^{\prime})\cos(\omega
t^{\prime}+\phi)e^{i\omega_{c}t^{\prime}}e^{-i\alpha_{0}(t-t^{\prime})}dt^{\prime}$
$\displaystyle=e^{-i\alpha_{0}t}\int_{-\infty}^{t}\Omega(t^{\prime})\cos(\omega
t^{\prime}+\phi)e^{i(\omega_{c}+\alpha_{0})t^{\prime}}dt^{\prime}.$ (13)

At this point, it is convenient to re-define $\theta$ to be a function of
$\omega_{c}$ as well as $t$ so that

$\theta(\omega_{c},t)=\int_{-\infty}^{t}\Omega(t^{\prime})\cos(\omega
t^{\prime}+\phi)e^{i\omega_{c}t^{\prime}}dt^{\prime}.$ (14)

With this notation, we can write (13) simply as

$\int_{-\infty}^{t}\dot{\theta}(t^{\prime})e^{\alpha(t^{\prime},t)}dt^{\prime}=e^{-i\alpha_{0}t}\theta(\omega_{c}+\alpha_{0},t).$
(15)

Substituting (15) back into (9), we have a closed-form solution for pulses of
arbitrary shapes:

$\tilde{f}_{1}(t)=-\frac{i}{2}[\theta(\omega_{c},t)+e^{-i\alpha_{0}t}\theta(\omega_{c}+\alpha_{0},t)],$
(16)

where $\alpha_{0}$ is a constant given by (12). From now on, tilde will be
used to denote the approximate solutions, while a notation without tilde
denotes an exact solution.

###  II.3 Sequence of Approximate Solutions

In order to obtain a more accurate solution of $\eqref{fnonlinear}$, we use
the recursive relation (8) to get a sequence of closed-form solutions
$\tilde{f}_{k}$ similar to (16). Taking the limit of $k\to\infty$ will lead us
to an improved solution.

First, note that the general solution to the linear equation (8) is:

$\displaystyle f_{k+1}$
$\displaystyle=e^{2i\int_{-\infty}^{t}f_{k}(t^{\prime})\dot{\theta}^{*}(t^{\prime})dt^{\prime}}\times$
(17)
$\displaystyle\left(\int_{0}^{t}e^{-2i\int_{-\infty}^{t^{\prime}}f_{k}(t^{\prime\prime})\dot{\theta}^{*}(t^{\prime\prime})dt^{\prime\prime}}[-if_{k}^{2}(t^{\prime})\dot{\theta}^{*}(t^{\prime})-i\dot{\theta}(t^{\prime})]dt^{\prime}\right),$

where the initial condition $f(0)=0$ has been used. Integrating the first
integrand by parts and bringing the outside exponential inside the integral
yield

$\displaystyle
f_{k+1}(t)=\frac{1}{2}\left(f_{k}(t)-\int_{0}^{t}e^{\beta_{k}(t^{\prime},t)}\left(\dot{f}_{k}(t^{\prime})+2i\dot{\theta}(t^{\prime})\right)dt^{\prime}\right),$
(18)

where

$\displaystyle\beta_{k}(t^{\prime},t)=e^{2i\int_{t^{\prime}}^{t}f_{k}(t^{\prime})\dot{\theta}^{*}(t^{\prime})dt^{\prime}}.$
(19)

Using the same idea as with (16), we approximate $\beta_{k}$ with a linear
function of $(t-t^{\prime})$

$\beta_{k}(t^{\prime},t)=-i\alpha_{k}(t-t^{\prime}),\;\;\alpha_{k}=-\frac{2}{\tau}\int_{0}^{\tau}\dot{\theta}^{*}(t^{\prime})\tilde{f}_{k}(t^{\prime})dt^{\prime}.$
(20)

Notice that we have used $\tilde{f}_{k}$, which is yet to be found, in the
expression for $\alpha_{k}$ above. Substituting this back into (18), and using
the same trick as in (15), we have

$\displaystyle\tilde{f}_{k+1}(t)=\frac{1}{2}\\{\tilde{f}_{k}(t)$
$\displaystyle-2ie^{-i\alpha_{k}t}\theta(\omega_{c}+\alpha_{k},t)$
$\displaystyle-\int_{0}^{t}e^{-i\alpha_{k}(t-t^{\prime})}\dot{\tilde{f}}_{k}(t^{\prime})dt^{\prime}\\},$
(21)

where we again note that all $f_{k}$ are tilded.

From (21), it is possible to derive an accurate solution for $\tilde{f}_{k+1}$
and then to take the limit $k\to\infty$. Although this would be the cleanest
way to derive the limiting solution, it is very tedious. Here, we offer a more
elegant approach. Assume that all functions in (21) are continuous and that
taking the limit commutes with integration. Take the limit $k\to\infty$ on
both sides and denote
$\tilde{f}_{\infty}=\displaystyle{\lim_{k\to\infty}}\tilde{f}_{k}$ and
$\alpha_{\infty}=\displaystyle{\lim_{k\to\infty}}\alpha_{k}$. We then have

$\displaystyle\frac{1}{2}\tilde{f}_{\infty}(t)=$ $\displaystyle-
ie^{-i\alpha_{\infty}t}\theta(\omega_{c}+\alpha_{\infty},t)$
$\displaystyle-\frac{1}{2}e^{-i\alpha_{\infty}t}\int_{0}^{t}e^{i\alpha_{\infty}t^{\prime}}\dot{\tilde{f}}_{\infty}(t^{\prime})dt^{\prime}.$
(22)

Introducing
$g(t)=e^{i\alpha_{\infty}t}\tilde{f}_{\infty}(t)+i\theta(\omega_{c}+\alpha_{\infty},t)$
and integrating the exponential by parts lead to

$\frac{1}{2}g(t)=-\frac{1}{2}\int_{0}^{t}\left(\dot{g}(t^{\prime})-i\alpha_{\infty}(g(t^{\prime})-i\theta(\omega_{c}+\alpha_{\infty},t^{\prime}))\right)dt^{\prime},$
(23)

and, after differentiation,

$\dot{g}(t)=\frac{i\alpha_{\infty}}{2}g(t)+\frac{\alpha_{\infty}}{2}\theta(\omega_{c}+\alpha_{\infty},t),$
(24)

which yields

$g(t)=e^{\frac{i\alpha_{\infty}}{2}t}\left(\int_{0}^{t}e^{\frac{-i\alpha_{\infty}}{2}t^{\prime}}\frac{\alpha_{\infty}}{2}\theta(\omega_{c}+\alpha_{\infty},t^{\prime})dt^{\prime}+C\right).$
(25)

Applying the initial condition $g(0)=0$ and integrating by parts, $g(t)$ is
simplified to

$g(t)=i\theta(\omega_{c}+\alpha_{\infty},t)-ie^{\frac{i\alpha_{\infty}}{2}t}\theta(\omega_{c}+\frac{\alpha_{\infty}}{2},t).$
(26)

Substituting $g(t)$ back to $\tilde{f}_{\infty}$, we have

$\tilde{f}_{\infty}(t)=-ie^{-i\lambda t}\theta(\omega_{c}+\lambda,t),$ (27)

where $\lambda=\alpha_{\infty}/2$ is a constant that we will find next.

The downfall of this shortcut derivation is that we do not immediately have an
expression for $\lambda$ (or $\alpha_{\infty}$). To find $\lambda$, let us
consider $\lambda$ as an arbitrary constant and try to identify the value of
$\lambda$ that would make (30) most accurate. Substituting (30) back into (5),
we find

$\lambda=i\dot{\theta}^{*}(\omega_{c},t)e^{-i\lambda
t}\theta(\omega_{c}+\lambda,t),$ (28)

which provides an equation the optimal $\lambda$ should satisfy. Notice that
the right hand side is time-dependent, so averaging it and using
$\dot{\theta}^{*}(\omega_{c},t)e^{i\nu t}=\dot{\theta}^{*}(\omega_{c}+\nu,t)$
lead us to the equation for the optimal $\lambda$,

$\lambda=\frac{i}{\tau}\int_{0}^{\tau}\dot{\theta}^{*}(\omega_{c}+\lambda,t^{\prime})\theta(\omega_{c}+\lambda,t^{\prime})dt^{\prime}.$
(29)

Admittedly, (29) does not provide a closed-form expression for $\lambda$.
However, the role of $\lambda$ is as a particularly suitable constant that
makes the approximate solution of the form (30) most accurate. The fact that
the form (30) is closed-form is the crucial point here.

Thus, we can summarize our solution for the Schrödinger equation as follows:

$\tilde{f}_{\infty}(t)=-ie^{-i\lambda t}\theta(\omega_{c}+\lambda,t),$ (30)

where

$\lambda=\frac{i}{\tau}\int_{0}^{\tau}\dot{\theta}^{*}(\omega_{c}+\lambda,t^{\prime})\theta(\omega_{c}+\lambda,t^{\prime})dt^{\prime},$
(31)

and

$\theta(\nu,t)=\int_{-\infty}^{t}\Omega(t^{\prime})\cos(\omega
t^{\prime}+\phi)e^{i\nu t^{\prime}}dt^{\prime}.$ (32)

##  III Gaussian Pulse Excitation

Figure 1: The percent $L^{2}$ error of the respective approximate solution as
compared to the exact numerical solution of (5). Top left: the zero-order
approximation (6). Top right: the approximate solution (16). Bottom left: the
limiting approximate solution (30). Bottom right: the exact solution of the
linear equation (7). The solid black line is the contour on which the error is
equal to 10%.

Let us study the conditions for which our solutions (16) and (30) are
applicable, as well as their accuracy.

First of all, it should be pointed out that the main result of Ref. [15] can
be obtained from our solution (16). For that, let us assume
$\alpha_{0}\ll\omega_{c}$ so that
$\theta(\omega_{c}+\alpha_{0},t)=\theta(\omega_{c},t)$. In this case, denoting
$\alpha_{0}=\eta^{2}\omega_{c}$, equation (16) becomes

$\tilde{f}_{1}(t)=-\frac{i}{2}\left(1+e^{-i\eta^{2}\omega_{c}t}\right)\theta(\omega_{c},t),$
(33)

which is exactly the equation (16) in Ref. [15]. The latter solution was shown
to be accurate for a square pulse under the condition [15]

$\left(\frac{\omega}{\omega_{c}}\right)^{2}+\left(\frac{\Omega_{0}}{\omega}\right)^{2}\ll
1.$ (34)

Next, we examine how accurate the solution (16) and its generalization (30)
are for more realistic pulses. As a quantity that expresses the accuracy of a
solution, we choose the $L^{2}$-norm of the deviation of the solution from the
exact solution. More specifically, the ratio of the latter quantity to the
$L^{2}$-norm of the exact solution itself - we call this ratio the relative
$L^{2}$ error. We numerically calculate this quantity for a Gaussian envelope
centered in the middle of the pulse

$\Omega(t)=\Omega_{0}\exp{-\frac{(t-\frac{\tau}{2})^{2}}{2\sigma^{2}}},$ (35)

where we mainly explore the dependence of the relative error described above
on the quantities $\Omega_{0}/\omega$ and $\omega_{c}/\omega$, since they were
shown in Ref. [15] to be the relevant parameters that determine applicability
of a solution.

In Fig. 1, we plot the relative error described above for the approximate
solutions (6), (16) and (30). As one can see, our most recent solution is the
most accurate of the three, especially in the region where the previous
solutions did not apply - for the high $\omega_{c}/\omega$ ratios. At the same
time, we can see that the applicability conditions are generally the same as
(34). Namely, the system must be far off resonance.

However, a slight modification of (34) is in order. Notice that our solution
applies not only for the situations when $\omega/\omega_{c}\ll 1$, but also
$\omega_{c}/\omega\gg 1$. As a result, the following condition is more
appropriate

$\min\left(\left(\frac{\omega_{c}}{\omega}\right)^{2},\left(\frac{\omega}{\omega_{c}}\right)^{2}\right)+\min\left(\left(\frac{\Omega_{0}}{\omega}\right)^{2},\left(\frac{\omega}{\Omega_{0}}\right)^{2}\right)\ll
1.$ (36)

As demonstrated in Fig. 1, our solution accurately predicts the behavior of
the TLS driven by a realistic Gaussian pulse. While the solution applies to
arbitrary pulse shapes, a Gaussian pulse offers an exclusively convenient
property under the current context, namely, the function
$\theta(\omega_{c},t)$ defined by (14) can be evaluated explicitly in closed-
form for a Gaussian pulse,

$\theta_{\text{Gaussian}}(\omega_{c},t)=\frac{\sigma\Omega_{0}}{4}\left(\zeta(t,\omega)e^{i\phi}+\zeta(t,-\omega)e^{-i\phi}\right),$
(37)
$\zeta(t,\omega)=e^{i\frac{\nu+\omega}{2}-\frac{1}{2}\sigma^{2}(\nu+\omega)^{2}}\left(\text{erf}\left(\frac{1+2i\nu\sigma^{2}+2i\sigma^{2}\omega}{2\sqrt{2}\sigma}\right)-\text{erf}\left(\frac{1-2t+2i\nu\sigma^{2}+2i\sigma^{2}\omega}{2\sqrt{2}\sigma}\right)\right).$
(38)

Such closed-form solutions have not been found for other common pulse shapes,
such as the Lorentzian or the hyperbolic secant pulses.

##  IV Discussion

Another insightful comparison of (30) is with the classical rotating-wave
approximation(RWA) regime. Under the RWA, the equations of motion (2) become

$\displaystyle\dot{C}(t)=-\frac{i}{2}\Omega(t)D(t),$ (39)
$\displaystyle\dot{D}(t)=-\frac{i}{2}\Omega(t)C(t).$ (40)

Introducing the area of the pulse
$A(t)=\int_{0}^{t}\Omega(t^{\prime})dt^{\prime}$, we can find a general
solution of the TLS in the form of

$\displaystyle C(t)=-i\sin\frac{A(t)}{2},$ (41) $\displaystyle
D(t)=\cos\frac{A(t)}{2}.$ (42)

Since $f(t)=C(t)/D(t)$, we have

$f=-i\tan\frac{A}{2}$ (43)

as the general solution under the RWA.

To appreciate the similarity of our solution, notice that we have the closed-
form solution (30) and the relation

$\theta(\omega_{c},t)=\frac{e^{i\phi}\tilde{A}(\omega_{c}+\omega,t)+e^{-i\phi}\tilde{A}(\omega_{c}-\omega,t)}{2}.$
(44)

Under RWA, $\tilde{A}(\omega_{c}+\omega,t)=0$ and
$\tilde{A}(\omega_{c}-\omega,t)=A(t)$, which means that up to a phase factor,
our solution gives

$\tilde{f}_{\infty}(t)=-i\frac{A}{2}.$ (45)

This suggests that, near resonance (which is outside of its region of
applicability), our solution corresponds to the linearization of the RWA
solution. This apparently is the result of multiple linearizations we have
taken in deriving $\tilde{f}_{k}$ to obtain the closed-form solution. With
this observation in mind, we aim to overcome the oversimplification of our
approximations and try to capture more of the nonlinear behaviors of the
original Schrödinger equation. One approach is outlined as follows.

Consider again the general equation (5) and use
$\dot{\theta}(t)=\Omega(t)\cos(\omega t+\phi)e^{i\omega_{c}t}$. Taking
$g(t)=\Omega(t)\cos(\omega t+\phi)$ and multiplying through by
$e^{-i\omega_{c}t}$, we arrive at

$e^{-i\omega_{c}t}\dot{f}(t)=ig(t)e^{-2i\omega_{c}t}f^{2}(t)-ig(t).$ (46)

Substituting $z(t)=e^{-i\omega_{c}t}f(t)$ into the above equation leads to

$\dot{z}(t)=ig(t)z^{2}(t)-i\omega_{c}z(t)-ig(t).$ (47)

Now expanding $z(t)$ as a perturbation series in $\omega_{c}$, we obtain for
the zeroth order

$\dot{z}_{0}(t)=ig(t)(z_{0}^{2}(t)-1),$ (48)

which gives

$z_{0}(t)=\cosh^{-1}\left(i\int_{0}^{t}g(t)dt\right).$ (49)

Meanwhile, the $n$th-order can be written as

$\dot{z}_{n}(t)=2ig(t)z_{n}(t)z_{0}(t)+i\displaystyle{\sum_{j=1}^{n-1}}z_{j}(t)z_{n-j}(t)-iz_{n-1}(t).$
(50)

The general solution to this linear equation, given the initial condition
$z(0)=0$, is:

$z_{n}(t)=ie^{w(t)}\int_{0}^{t}e^{-w(t^{\prime})}\left(\displaystyle{\sum_{j=1}^{n-1}}z_{j}(t^{\prime})z_{n-j}(t^{\prime})-z_{n-1}(t^{\prime})\right)dt^{\prime},$
(51)

where $w(t)=2i\int_{0}^{t}g(t^{\prime})z_{0}(t^{\prime})dt^{\prime}$. Together
with (49), (51) defines a sequence of functions that converge to the solution,
similar to the sequence $f_{k}(t)$ discussed in Section II.

Note that, compared to the linearized sequence $f_{k}(t)$, the sequence of
$z_{n}(t)$ is a much closer approximation to the accurate solution. Even the
first-order approximation (49) captures the nonlinear behavior of the system.
If a sequence of closed-form approximations to (51) could be found, similar to
the sequence $\tilde{f}_{k}$ presented in Section II, the limit of such
sequence $z_{\infty}$ would be a promising approximation for capturing the
nonlinear behaviors of the system. It seems plausible that a combination of
such solution with the solution (30) would allow us to extend the conditions
of applicability of (30) to the case when
$\omega\sim\omega_{c},\;\omega\sim\Omega_{0}$, and capture RWA better than a
simple linear approximation. However, so far, we have been unable to find any
approach to obtain a closed-form solution in this way. The difficulty is
precisely the nonlinear features embedded into each (51).

##  V Conclusion

In this work, we have developed a closed-form analytical solution of the
Schrödinger equation that describes a TLS driven by a far-off-resonance few-
cycle pulse without using the RWA. Unlike our previous solution, which only
holds for square pulses [15], the current solution is generalized to arbitrary
pulse shapes. We have identified the conditions under which our solution
accurately predicts the behaviors of the system, and have demonstrated that
the new solution offers improved accuracy compared to the previous solution
under similar conditions. We have also applied the general solution to
Gaussian pulses as a demonstration of its applicability. Finally, we suggest a
possible alternative approach that can lead to a more accurate solution by
capturing the nonlinear behaviors of the system and extending the
applicability of the solution to the case of
$\omega\sim\omega_{c},\;\omega\sim\Omega_{0}$. This work lays out a potential
pathway toward an analytical theory for ultrafast QCC beyond the RWA.

## References

  * Brif  _et al._ [2010] C. Brif, R. Chakrabarti, and H. Rabitz, Control of quantum phenomena: past, present and future, New Journal of Physics 12, 075008 (2010).
  * Genkin [1998] G. M. Genkin, Rabi frequency and nonlinearity of a two-level atom for an ultrashort optical pulse, Phys. Rev. A 58, 758 (1998).
  * Jirauschek  _et al._ [2005] C. Jirauschek, L. Duan, O. D. Mücke, F. X. Kärtner, M. Wegener, and U. Morgner, Carrier-envelope phase-sensitive inversion in two-level systems, J. Opt. Soc. Am. B 22, 2065 (2005).
  * Xie  _et al._ [2009] X.-T. Xie, M. Macovei, M. Kiffner, and C. H. Keitel, Probing quantum superposition states with few-cycle laser pulses, J. Opt. Soc. Am. B 26, 1912 (2009).
  * Kumar and Sarma [2012] P. Kumar and A. K. Sarma, Gaussian and sinc-shaped few-cycle-pulse-driven ultrafast coherent population transfer in $\Lambda$-like atomic systems, Phys. Rev. A 85, 043417 (2012).
  * Došlić [2006] N. Došlić, Generalization of the rabi population inversion dynamics in the sub-one-cycle pulse limit, Phys. Rev. A 74, 013402 (2006).
  * Roudnev and Esry [2007] V. Roudnev and B. D. Esry, General theory of carrier-envelope phase effects, Phys. Rev. Lett. 99, 220406 (2007).
  * Wu and Yang [2007] Y. Wu and X. Yang, Strong-coupling theory of periodically driven two-level systems, Phys. Rev. Lett. 98, 013601 (2007).
  * Rostovtsev  _et al._ [2009] Y. V. Rostovtsev, H. Eleuch, A. Svidzinsky, H. Li, V. Sautenkov, and M. O. Scully, Excitation of atomic coherence using off-resonant strong laser pulses, Phys. Rev. A 79, 063833 (2009).
  * Jha and Rostovtsev [2010a] P. K. Jha and Y. V. Rostovtsev, Analytical solutions for a two-level system driven by a class of chirped pulses, Phys. Rev. A 82, 015801 (2010a).
  * Jha  _et al._ [2010] P. K. Jha, H. Eleuch, and Y. V. Rostovtsev, Coherent control of atomic excitation using off-resonant strong few-cycle pulses, Phys. Rev. A 82, 045805 (2010).
  * Jha and Rostovtsev [2010b] P. K. Jha and Y. V. Rostovtsev, Coherent excitation of a two-level atom driven by a far-off-resonant classical field: Analytical solutions, Phys. Rev. A 81, 033827 (2010b).
  * Beg [2020] Magnus expansion applied to a dissipative driven two-level system, Results in Physics 17, 103098 (2020).
  * Arkhipov  _et al._ [2019] R. M. Arkhipov, A. V. Pakhomov, M. V. Arkhipov, I. Babushkin, A. Demircan, U. Morgner, and N. N. Rosanov, Unipolar subcycle pulse-driven nonresonant excitation of quantum systems, Opt. Lett. 44, 1202 (2019).
  * Zeng and Duan [2021] B. Zeng and L. Duan, An analytical theory of cep-dependent coherence driven by few-cycle pulses, Physics Letters A 398, 127292 (2021).
  * Solá  _et al._ [1999] I. R. Solá, V. S. Malinovsky, B. Y. Chang, J. Santamaria, and K. Bergmann, Coherent population transfer in three-level $\Lambda$ systems by chirped laser pulses: Minimization of the intermediate-level population, Phys. Rev. A 59, 4494 (1999).
  * Cai  _et al._ [2013] X. Cai, J. Zhao, Z. Wang, and Q. Lin, Ultrafast coherent population transfer in two- and three-level quantum systems using sub-cycle and single-cycle pulses, Journal of Physics B Atomic Molecular Physics 46, 175602 (2013).
  * Scully  _et al._ [2008] M. O. Scully, Y. V. Rostovtsev, A. Svidzinsky, and J. Chang, Xuv coherent raman superradiance, Journal of Modern Optics 55, 3219 (2008).
  * Leb [2013] Models of few optical cycle solitons beyond the slowly varying envelope approximation, Physics Reports 523, 61 (2013), models of few optical cycle solitons beyond the slowly varying envelope approximation.
  * Hughes [1998] S. Hughes, Breakdown of the area theorem: Carrier-wave rabi flopping of femtosecond optical pulses, Phys. Rev. Lett. 81, 3363 (1998).
  * Novitsky [2012] D. V. Novitsky, Propagation of subcycle pulses in a two-level medium: Area-theorem breakdown and pulse shape, Phys. Rev. A 86, 063835 (2012).
  * Brown and Meath [1998] A. Brown and W. J. Meath, On the effects of absolute laser phase on the interaction of a pulsed laser with polar versus nonpolar molecules, The Journal of Chemical Physics 109, 9351 (1998), https://doi.org/10.1063/1.477596 .
  * Kumar and Sarma [2013] P. Kumar and A. K. Sarma, Frequency-modulated few-cycle optical-pulse-train-induced controllable ultrafast coherent population oscillations in two-level atomic systems, Phys. Rev. A 87, 025401 (2013).
  * Yang  _et al._ [2014] W.-X. Yang, A.-X. Chen, H. Guo, Y. Bai, and R.-K. Lee, Carrier-envelope phase control electron transport in an asymmetric double quantum dot irradiated by a few-cycle pulse, Optics Communications 328, 96 (2014).
  * Jha  _et al._ [2014] P. K. Jha, H. Eleuch, and F. Grazioso, Ultra-short strong excitation of two-level systems, Optics Communications 331, 198 (2014).


# Static and dynamic magnetic properties of the spin-$\frac{5}{2}$ triangle
lattice antiferromagnet Na3Fe(PO4)2 studied by 31P NMR

Devi V. Ambika1,2, Qing-Ping Ding1,2, Sebin J. Sebastian3, Ramesh Nath3, and
Yuji Furukawa1,2 1 Ames National Laboratory, U.S. DOE, Iowa State University,
Ames, Iowa 50011, USA  2 Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  3 School of Physics, Indian Institute of
Science Education and Research Thiruvananthapuram-695551, India
<EMAIL_ADDRESS>
###### Abstract

31P nuclear magnetic resonance (NMR) measurements have been carried out to
investigate the magnetic properties and spin dynamics of Fe3+ ($S$ = 5/2)
spins in the two-dimensional triangular lattice (TL) compound Na3Fe(PO4)2. The
temperature ($T$) dependence of nuclear spin-lattice relaxation rates
($1/T_{1}$) shows a clear peak around Néel temperature, $T_{\rm N}=10.9$ K,
corresponding to an antiferromagnetic (AFM) transition. From the temparature
dependence of NMR shift ($K$) above $T_{\rm N}$, an exchange coupling between
Fe3+ spins was estimated to be $J/k_{\rm B}\simeq 1.9$ K using the spin-5/2
Heisenberg isotropic-TL model. The temperature dependence of $1/T_{1}T$
divided by the magnetic susceptibility ($\chi$), $1/T_{1}T\chi$, above $T_{\rm
N}$ proves the AFM nature of spin fluctuations below $\sim$ 50 K in the
paramagnetic state. In the magnetically ordered state below $T_{\rm N}$, the
characteristic rectangular shape of the NMR spectra is observed, indicative of
a commensurate AFM state in its ground state. The strong temperature
dependence of 1/$T_{1}$ in the AFM state is well explained by the two-magnon
(Raman) process of the spin waves in a 3D antiferromagnet with a spin-
anisotropy energy gap of 5.7 K. The temperature dependence of sublattice
magnetization is also well reproduced by the spin waves. Those results
indicate that the magnetically ordered state of Na3Fe(PO4)2 is a conventional
3D AFM state, and no obvious spin frustration effects were detected in its
ground state.

##  1 Introduction

Geometrically frustrated magnets have attracted great interest because the
systems show peculiar properties such as the suppression of the magnetic
ordering, unconventional spin states like the spin liquid [1], and large
fluctuations in the ground state. Initially, the ground state of two-
dimensional (2D) triangular lattice (TL) antiferromagnets was predicted to
host a resonating valence bond state [2], and the long-range magnetic ordering
in TL was identified to occur when the spins orient at 120∘ angle relative to
neighboring spins [3, 4]. Later studies of novel compounds proved the
existence of non-trivial ground states in such systems. For example, the
presence of a weak easy-axis anisotropy in $S\geq 3/2$ triangular lattice
antiferromagnets leads to a collinear ordered state [5].

Figure 1:  (a) Crystal structure of Na3Fe(PO${}_{4})_{2}$ showing corner-
shared FeO6 octahedra and PO4 tetrahedra forming triangular layers. Na atoms
are located in-between the layers. (b) Fe3+ ions connected via PO4 tetrahedra
forming the triangular lattice in the $ab$ plane.

Na3Fe(PO${}_{4})_{2}$ is one of the series compounds of the family
$AA^{\prime}M$(RO${}_{4})_{2}$, (A = Ba, K, Rb, Na; $A^{\prime}$ =
$\rm{Na}_{2}$, $\rm{Ag}_{2}$; M = Cr, Mn, Fe, Co, Ni, and Cu; and R = P, V,
As) [6, 7, 8] which provides a good platform to study systematically the
magnetic properties of triangular lattice antiferromagnets with $S\geq 3/2$.
Depending upon the size of the non-magnetic spacers and exchange coupling, the
systems can host either a magnetically ordered phase or a spin liquid phase.
Na3Fe(PO${}_{4})_{2}$ crystallizes in a monoclinic space group ($C2/c$) with
the lattice parameters of $a$ = 9.0714(2) Å, $b$ = 5.032(1) Å, and $c$ =
13.8683(3) Å at room temperature [9, 10]. Fe3+ ions ($S$ = 5/2) in distorted
FeO6 octahedra form an $S$ = 5/2 anisotropic TL via corner-sharing with PO4
tetrahedra in the $ab$ plane as shown in Fig. 1, where the triangles are
slightly distorted due to the two different Fe-Fe distances of 5.182 Å and
5.028 Å [11]. Recently, Sebastian $et~{}al.$ reported that
Na3Fe(PO${}_{4})_{2}$ exhibits an antiferromagnetic (AFM) ordering with a Néel
temperature of $T_{\rm N}=10.6$ K with a reduced Fe ordered magnetic moment of
$\sim$ 1.52 $\mu_{\rm B}$ at 1.6 K in zero magnetic field [11]. The average
exchange coupling between Fe3+ spins was estimated to be $J/k_{\rm B}\simeq
1.8$ K by fitting magnetic susceptibility $\chi(T)$ data using the high-
temperature series expansion for a spin-5/2 Heisenberg TL model. The authors
also determined, from neutron diffraction measurements, a collinear AFM
structure with a propagation vector $q$ = (1, 0, 0) where the Fe spins are
coupled ferromagnetically along the $b$ axis or in the [100] plane [11].
However, the detailed magnetic properties, especially magnetic fluctuations,
of Na3Fe(PO${}_{4})_{2}$ have still not been investigated.

Nuclear magnetic resonance (NMR) is a very powerful experimental tool to probe
the static and dynamic magnetic properties of materials from a microscopic
point of view. In this paper, we report the results of 31P NMR measurement on
the polycrystalline Na3Fe(PO${}_{4})_{2}$ sample in the temperature range of
$T=1.5-300$ K. The NMR spectra, nuclear spin-lattice relaxation rate 1/$T_{1}$
and nuclear spin-spin relaxation rate 1/$T_{2}$ measurements show the AFM
long-range order below $T_{\rm N}\sim$ 11 K. In the paramagnetic state above
$T_{\rm N}$, from the analysis of NMR shift $K(T)$ using the high-temperature
series expansion, $J/k_{\rm B}$ was estimated to be $\simeq 1.9$ K, which is
consistent with the previously reported value of $J/k_{\rm B}\simeq 1.8$ K.
From the temperature dependences of 1/$T_{1}$ and 1/$T_{2}$ in the
paramagnetic state above $T_{\rm N}$, AFM spin fluctuations were found to be
developed below $\sim$ 50 K. In the AFM ordered state, we observed a
rectangular shape of the spectrum, evidencing a commensurate AFM state.
However, additional peaks in the spectra in the AFM state were observed above
3 T, which are probably due to the spin flop reported previously [11]. The
temperature dependences of the sublattice magnetization and 1/$T_{1}$ in the
AFM state were well explained by spin waves for a conventional three-
dimensional antiferromagnet with a finite magnetic anisotropy gap $\sim$ 5.7
K, indicating that the AFM state is rigid without any obvious effects of spin
frustrations in its ground state.

##  2 Experimental

Polycrystalline Na3Fe(PO${}_{4})_{2}$ sample was synthesized using the
conventional solid-state method [11]. NMR measurements of 31P ($I$ =
$\frac{1}{2}$, $\frac{\gamma_{N}}{2\pi}$ = $17.237$ MHz/T) nuclei were
conducted using a laboratory-built phase-coherent spin-echo pulse
spectrometer. The 31P NMR spectra were obtained by sweeping magnetic field $H$
at fixed resonance frequencies. 31P zero-shift position ($H_{\rm ref}$) for
each resonance frequency was determined by the 31P NMR of the nonmagnetic
reference sample H3PO4. The 31P spin-lattice relaxation rates $1/T_{1}$ were
measured by the saturation recovery method. $1/T_{1}$ at each temperature was
determined by fitting the nuclear magnetization $M(t)$ versus time $t$ using
the single exponential function

$1-\frac{M(t)}{M(\infty)}=A~{}exp\left[-\left(\frac{t}{T_{1}}\right)\right],$
(1)

where $M(\infty)$ is the equilibrium nuclear magnetization at
$t\rightarrow\infty$.

Nuclear spin-spin relaxation times $T_{2}$ were determined by fitting spin-
echo decay curves using the following equation:

$M(2\tau)=M_{0}~{}exp\left[-\left(\frac{2\tau}{T_{2}}\right)^{\beta}\right],$
(2)

where $\tau$ is the time between $\pi$/2 and $\pi$ pulses. In the paramagnetic
state, Lorentzian relaxation behavior with $\beta=1$ was observed. However, in
the AFM state, $\beta$ increases with decreasing temperature and becomes
$\approx 2$ at low temperatures below $T=3.4$ K, showing Gaussian decay
behavior.

##  3 Results and discussions

###  3.1 31P NMR spectrum in the paramagnetic state

Figure 2:  Field-swept 31P NMR spectra at different temperatures measured at
$f$ = 32.5 MHz. The vertical dash-dotted line corresponds to the 31P zero-
shift or non-magnetic reference field ($H_{\rm ref}\approx 1.886$ T). The
inset shows the spectrum at 40 K with the simulated spectrum (red curve) with
$K_{\rm iso}\simeq 4.91$%, $K_{\rm ax}\simeq-0.28$%. Figure 3: Temperature-
dependent 31P NMR iso-shift $K_{\rm iso}$ and $K_{\rm ax}$ measured at 32.5
MHz. The red curve shows the result of the fit to $K_{\rm iso}$ using Eq. 3,
above 20 K. Inset: $K$ vs. $\chi$ plots for $K_{\rm iso}$ and $K_{\rm ax}$
where we used $\chi$ data reported in Ref. [11]. The thick solid lines are
fitting results.

Figure 2 shows the temperature dependence of the field-swept 31P NMR spectra
measured on the polycrystalline sample at a frequency of f = 32.5 MHz ($H_{\rm
ref}\approx 1.886$ T) at several temperatures ($T>T_{\rm N}$). The observed
spectra show a characteristic shape of the anisotropic powder pattern. With
decreasing temperature, the spectra become broader and shift to lower magnetic
fields. Then the spectra move to slightly higher magnetic fields below 20 K.

The asymmetric shape of spectra was well reproduced by the calculated spectra
where we introduced an anisotropy in NMR shift $K$. The red curve in the inset
of Fig. 2 shows the typical calculated powder-pattern spectrum with isotropic
and axial components of NMR shifts ($K_{\rm iso}$ and $K_{\rm ax}$) and
reproduced the observed spectrum (open black circles) very well. Here the NMR
shift $K$ is described by $K$ = $K_{\rm iso}$ \+ $K_{\rm ax}$(3cos2$\theta$
-1) where $\theta$ is the angle between the principal axis of hyperfine tensor
at the P site and the external magnetic field [12].

From the fitting of the spectra measured at different temperatures, we
determined the temperature dependences of $K_{\rm iso}$ and $K_{\rm ax}$ as
shown in Fig. 3. The temperature dependences of $K_{\rm iso}$ and $K_{\rm ax}$
are essentially the same physically, although their sign and magnitude are
different. The absolute values of $K_{\rm iso}$ and $K_{\rm ax}$ increase with
decreasing temperature in the high-temperature region following the Curie-
Weiss behavior. Further decrease in temperature leads to a broad maximum
(minimum) in $K_{\rm iso}$ ($K_{\rm ax}$) around 17 K, analogous to the
$\chi({T})$ data depicting a short-range AFM ordering. $K$ has a direct
relation with the spin susceptibility $\chi_{\rm spin}$, and is given by the
equation

$K=K_{0}+\frac{A_{\rm hf}}{N_{\rm A}\mu_{\rm B}}\chi_{\rm spin},$ (3)

where $K_{0}$ is the temperature-independent chemical shift, $N_{A}$ is the
Avogadro number, and $A_{\rm hf}$ is the hyperfine coupling constant between
the 31P nuclei and the Fe3+ electronic spins. Thus, the slope from the linear
fit of $K$ vs. $\chi$ plot with temperature as an implicit parameter yields
the corresponding $A_{\rm hf}$. As shown in the inset of Fig. 3, the linear
dependence of $K_{\rm iso}$ and $K_{\rm ax}$ on the $\chi$ over the whole
temperature range ($T>T_{\rm N}$) can be seen where the $\chi$ data are
obtained from Ref. [11], showing again the Curies-Weiss behavior in $K_{\rm
iso}$ and $K_{\rm ax}$ in the temperature region. The slopes yield the
isotropic part of the hyperfine coupling constant $A^{\rm iso}_{\rm
hf}~{}\simeq(4.29\pm 0.02)$ kOe/$\mu_{\rm B}$ and the axial part of the
hyperfine coupling constant $A^{\rm ax}_{\rm hf}~{}\simeq(-0.29\pm 0.01)$
kOe/$\mu_{\rm B}$.

The isotropic part of the $A_{\rm hf}$ originates from the transferred
hyperfine field due to the P($3s$)-O($2p$)-Fe(3$d$) covalent bond. On the
other hand, the axial part of the $A_{\rm hf}$ could be due to an anisotropic
hyperfine field at P atoms though the P($3p$)-O($2p$)-Fe(3$d$) covalent bond
and/or classical dipolar field from Fe3+ ions. From our calculation for the
classical dipolar field, it is found that the axial component of the classical
dipolar field is $|$0.18$|$ kOe/$\mu_{\rm B}$, the same order as the observed
$A^{\rm ax}_{\rm hf}~{}\simeq(-0.29\pm 0.01)$ kOe/$\mu_{\rm B}$, indicating
that the $A^{\rm ax}_{\rm hf}$ mainly comes from the classical dipolar field.

To estimate the exchange coupling $J$ between Fe3+ spins, we fitted the
temperature dependence of $K_{\rm iso}$ above 20 K by Eq. 3, where we used the
high-temperature series-expansion of $\chi_{\rm spin}$ for a spin-$S$
Heisenberg isotropic triangle lattice model [13] given by,

$\frac{N_{A}\mu_{B}^{2}g^{2}}{3\left|J\right|\chi_{\rm
spin}}=x+4+\frac{3.20}{x}-\frac{2.186}{x^{2}}+\frac{0.08}{x^{3}}+\frac{3.45}{x^{4}}-\frac{3.99}{x^{5}},$
(4)

with $x=k_{\rm B}T/|J|S(S+1)$. The Eq. 4 predicts the $\chi_{\rm spin}$ of a
triangular lattice accurately for $T\geq(J/k_{\rm B})S(S+1)$ [14]. With using
$A^{\rm iso}_{\rm hf}~{}\simeq 4.29$ kOe/$\mu_{\rm B}$, we fitted the data as
shown by the red curve in Fig. 3 and estimated $J/k_{\rm B}\simeq 1.9$ K, and
Landè g-factor $g~{}\simeq 2.18$. The exchange coupling obtained by this fit
is close to the previously reported value of $J/k_{\rm B}$ = 1.8 K from
$\chi(T)$ analysis [11].

###  3.2 31P NMR spectrum in the AFM state

Figure 4: Temperature-dependent 31P NMR spectra measured at 32.5 MHz, below
$T_{\rm N}$. The $H_{\rm int}$ broadens the spectra below $T_{\rm N}$, taking
a rectangular shape. The red curves show the calculated spectra using Eq. 6
for a commensurate AFM ordered state.

The 31P spectra measured at 32.5 MHz in the AFM state are shown in Fig. 4. The
spectra broaden below $T_{\rm N}$ and exhibit nearly rectangular shapes at low
temperatures. The immediate line broadening below $T_{\rm N}$ is due to the
static internal field at the P site produced by the ordered moments of Fe3+
ions in the magnetically ordered state.

In general, rectangular-shape NMR spectra are expected in rigid
antiferromagnets in powder samples where the angle between the internal field
$H_{\rm int}$ and the external magnetic field is randomly distributed [15].
These rectangular spectra were observed in several other compounds showing a
commensurate AFM ordered state [16].

Assuming a uniform $H_{\rm int}$ in a commensurate AFM ordered state, one can
express the NMR spectra, $f(H)$, as [17, 18]

$f(H)\propto\frac{H^{2}-H^{2}_{\rm int}+\omega^{2}_{\rm N}/\gamma^{2}_{\rm
N}}{H_{\rm int}H^{2}},$ (5)

where $\omega_{\rm N}/\gamma_{\rm N}=H_{0}$ is the central field, $\omega_{\rm
N}$ is the NMR frequency, $\gamma_{\rm N}$ is the nuclear gyromagnetic ratio,
and $H_{\rm int}$ is the internal field. The cut-off fields are represented as
$H_{0}+H_{\rm int}$ and $H_{0}-H_{\rm int}$, which produces the two sharp
edges of the spectrum. Due to the small distribution of internal fields, the
internal field $H_{\rm int}$ is modeled by a Gaussian distribution function to
smoothen the sharp edges. To obtain the spectra, the function $f(H)$ and the
distribution function $g(H)$ were convoluted as

$F(H)=\int_{0}^{\infty}f(H-H^{\prime})g(H^{\prime})\,dH^{\prime},$ (6)

where the Gaussian distribution function is given by

$g(H)=\frac{1}{\sqrt{2\pi\Delta H^{2}_{\rm int}}}exp\left(-\frac{(H-H_{\rm
int})^{2}}{2\Delta H^{2}_{\rm int}}\right).$ (7)

The red curves in Fig. 4 show calculated spectra, which reproduce the observed
spectra well. This indicates a commensurate AFM magnetic order below $T_{\rm
N}$, consistent with neutron diffraction measurements [11]. It is noted that
there is a small deviation between the calculated and observed spectra around
the edge of the spectrum at the high magnetic field side. It is probably due
to the additional signal intensity due to 23Na NMR signals (not shown) and/or
it may be due to partial orientation of the powder samples originating from
the application of magnetic field.

Figure 5:  Temperature dependence of the internal field $H_{\rm int}$ obtained
from the NMR spectra in the AFM state. The blue triangles and black closed
circles show the $H_{\rm int}$ estimated from fittings and from the half of
full-width at 25 % of the peak intensity of spectra, respectively. The green
solid curve is the fit of $H_{\rm int}$ vs. $T$ using Eq. 8 in the critical
regime, and the green dotted line exhibits the constant term of Eq. 8 in the
paramagnetic state. The red solid curve shows the calculated result using Eq.
12 for the magnon model in the AFM state.

Figure 5 shows the temperature dependence of $H_{\rm int}$ determined by the
aforementioned fittings of the spectra (blue triangles). It was difficult to
precisely determine the $H_{\rm int}$ near $T_{\rm N}$ by the fittings because
the commensurate AFM state is not well developed. Therefore we also estimated
the $H_{\rm int}$ by taking half of full-width at $25\%$ of the peak intensity
of the spectra shown by the black closed circles in order to determine $H_{\rm
int}$ near $T_{\rm N}$. As seen, both methods produce almost the same $H_{\rm
int}$ at $T$ $<<$ $T_{\rm N}$. As we approach the $T_{\rm N}$, $H_{\rm int}$
decreases sharply. Since $H_{\rm int}$ is proportional to the Fe3+ sublattice
magnetization, $H_{\rm int}(T)$ close to $T_{\rm N}$ is related to the
critical exponent ($\beta$) of the order parameter. To estimate the value of
$\beta$, we fit the $H_{\rm int}$ obtained by the line width in the critical
regime, very close to $T_{\rm N}$ using the following equation

$H_{\rm int}(T)=H_{\rm int}(0)\left(1-\frac{T}{T_{N}}\right)^{\beta}+c,$ (8)

where the constant value of $c$ is included to consider the line width in the
paramagnetic state. Our fit in the temperature range 10-10.6 K by Eq. 8 yields
$\beta\simeq 0.29$ with $H_{\rm int}(0)\simeq 1.1$ T, $T_{\rm N}=10.6$ K and
$c\simeq 0.013$ T. The obtained value of $\beta$ is slightly greater than the
expected value for 2D XY model ($\beta\simeq 0.231$) [19] and 2D Ising model
($\beta\simeq 0.25$) [20, 21] and is close to or is slightly smaller than the
expected values for 3D system (0.31-0.326 for 3D Ising model, 0.31-0.345 for
3D XY model and 0.33-0.367 for 3D Heisenberg model [22]). On the other hand,
when we change the fitting range to 9-10.6 K from 10-10.6 K, the $\beta$ value
decreases to 0.23 which is close to the value for the 2D XY model. Since
$\beta$ = 0.29 is obtained in the temperature regime close to $T_{\rm N}$, our
results may suggest that the spin dimensionality of the system very close to
$T_{\rm N}$ is most likely characterized by the 3D nature in NaFe3(PO4)2.

Figure 6:  31P NMR spectra measured at $T=4.3$ K for different applied
magnetic fields $H$, above and below the spin-flop transition. The arrows
depict the evolution of extra features with the field above the spin-flop
transition.

Finally, we show the field-dependent 31P NMR spectra measured at $T=4.3$ K in
Fig. 6. Interestingly additional peaks were observed with the increase in $H$
above $\sim 3$ T, whose positions seem to move to the center of the spectra
with increasing $H$. These results indicate the appearance of the P site with
different $H_{\rm int}$ under magnetic fields. It is interesting to point out
that a spin flop has been reported around 3.2 T at 4.2 K from the
magnetization measurements [11]. Although the origin of those additional peaks
in the NMR spectrum is not clear at present, if it were due to the spin flop,
our experimental results suggest that the spin flop may take place partially,
not uniformly, since the observed spectra consist of the rectangular powder
pattern due to the uniform internal field and the additional peaks with
different internal fields. Further experiments such as neutron diffraction
measurements under magnetic field and/or NMR measurements on single crystals
are required to clarify this point.

###  3.3 31P spin-lattice relaxation rate in the paramagnetic state

Figure 7:  Temperature dependence of $1/T_{1}$ measured at three frequencies
(32.5, 78.5, and 119 MHz). The red solid line represents the fits using Eq. 10
based on the model of the two magnon Raman process with $C$ = 2.5 (1/sK3) and
$T_{\rm AE}$ = 5.7 K in the AFM state. The green dotted line shows the
calculated result with the model for triangular antiferromagnet using Eq. 11
with $C^{\prime}$ = 1.3$\times$10-3 (1/sK5) and $T_{\rm AE}$ = 5.7 K. Inset
shows the $H$ dependence of $1/T_{1}$ at several temperatures (20, 40, 100,
and 300 K). The solid lines except for 20 K data are the expected $H$
dependence of 1/$T_{1}$ based on the 1D spin-diffusion model. It is noted that
1/$T_{1}$ data at 20 K shows a nearly $H$ independent behavior, indicating no
obvious spin diffusion effects. The green dotted curve shows the log(1/$H$)
dependence of 1/$T_{1}$ expected for the 2D spin-diffusion model fitted the
data at 300 K.

The spin-lattice relaxation rate $1/T_{1}$ for 31P nuclei was measured at the
peak position of the spectrum to understand the dynamical properties of Fe3+
spins. Figure 7 shows the temperature dependence of $1/T_{1}$ measured at $f$
= 32.5 MHz. Above $\sim$50 K, $1/T_{1}$ is nearly independent of temperature
and can be explained by the random fluctuation of paramagnetic moments in the
limit $T\gg J/k_{\rm B}$ [23]. As the temperature is lowered below $T$ $<$
$\sim$50 K, $1/T_{1}$ slowly decreases, showing a local minimum around 15 K,
and then peaks at $T_{\rm N}\simeq 10.9$ K, due to the magnetic long-range
order. Below $T_{\rm N}$, the $1/T_{1}$ decreases strongly with decreasing
temperature.

In the paramagnetic state, 1/$T_{1}$ depends on the resonance frequency
(magnetic field), as shown in Fig. 7. 1/$T_{1}$ decreases with an increase in
frequency in the high-temperature region. The difference in 1/$T_{1}$ is
narrowed down as the temperature is lowered, and below $T\sim$ 20 K, the data
sets in different frequencies nearly overlap with each other. It is known that
the long-wavelength ($q$ = 0) spin fluctuations in a Heisenberg magnet often
show diffusive dynamics. This is due to the divergence behavior of the
spectral density of the spin-spin correlations as $\omega$ goes to zero, which
depends on the dimensionality of the spin system. In 1D spin chains, such a
spin diffusion leads to a $1/\sqrt{H}$ field dependence of 1/$T_{1}$ as
observed in Sr2CuO3 [24] and (CH3)4NMnCl3 [25]. In the case of 2D materials,
$1/T_{1}$ is known to vary as log$(1/H)$ [26, 27, 28]. In the inset of Fig. 7,
we plot $1/T_{1}$ against $H$ at several temperatures. We fitted the
experimental data at 300 K with the relations 1/$T_{1}$ $\propto$ $1/\sqrt{H}$
and 1/$T_{1}$ $\propto$ $\log{1/H}$ which are shown by the green solid and
dotted curves in the inset of Fig. 7. As both $H$ dependences of $1/T_{1}$
capture the experimental observations, we do not conclude whether the system
is close to the 1D or 2D system within our experimental uncertainty. For other
temperatures such as 40 K and 100 K, we tentatively show the fits for the 1D
model. It is noted that 1/$T_{1}$ data at 20 K shows a nearly $H$-independent
behavior, indicating no obvious spin diffusion effects. In any case, we
consider that the $H$-dependent behavior of 1/$T_{1}$ observed at high
temperatures in NaFe3(PO4)2 originates from the spin diffusion effect, a
characteristic of low-dimensional spin systems.

Figure 8:  Temperature dependence of $1/T_{1}T\chi$ measured at $f$ = 32.5
MHz.

Finally, to see spin fluctuation effects in the paramagnetic state above
$T_{\rm N}$, it is useful to replot the data by changing the vertical axis
from 1/$T_{1}$ to $1/T_{1}T\chi$ as shown in Fig. 8. $1/T_{1}T$ is related to
the imaginary part of dynamic susceptibility
$\chi^{\prime\prime}\left(q,\omega_{\rm n}\right)$ at the NMR frequency
$\omega_{\rm N}$ as [29]

$\frac{1}{T_{1}T}=\frac{2\gamma_{N}^{2}k_{\mathrm{B}}}{N_{\mathrm{A}}^{2}}\sum_{q}|A(q)|^{2}\frac{\chi^{\prime\prime}\left(q,\omega_{\rm
N}\right)}{\omega_{\rm N}},$ (9)

where $k_{\mathrm{B}}$ is the Boltzmann’s constant. In Eq. 9, the sum is over
the wave vectors $\vec{q}$ within the first Brillouin zone, and $A(\vec{q})$
is the form factor of the hyperfine interactions. The uniform static
susceptibility $\chi=\chi^{\prime}(0,0)$, is the real part
$\chi^{\prime}\left(q,\omega_{\rm N}\right)$ with $q=0$ and $\omega_{\rm
N}=0$. Thus, by plotting $1/T_{1}T\chi$ vs. $T$, we get the temperature
dependence of $\sum_{q}|A(q)|^{2}\chi^{\prime\prime}\left(q,\omega_{\rm
N}\right)$ compared to that of $\chi^{\prime}(0,0)$, the uniform
susceptibility.

As shown in Fig. 8, $1/T_{1}T\chi$ is nearly independent of temperature in the
high temperature region down to $T\sim 50$ K. Thus, in this temperature
region, $\sum_{q}|A(q)|^{2}\chi^{\prime\prime}\left(q,\omega_{\rm N}\right)$
scales to $\chi^{\prime}(0,0)$. Below $T\sim$ 50 K till $T_{\rm N}$, it
increases with decreasing temperature, implying that
$\sum_{q}|A(q)|^{2}\chi^{\prime\prime}\left(q,\omega_{\rm N}\right)$ increases
more than $\chi^{\prime}(0,0)$. This result evidences a growth of spin
fluctuations with $q\neq 0$, attributing to AFM correlations.

###  3.4 31P spin-lattice relaxation rate in the AFM state

We observe a strong temperature dependence of spin-lattice relaxation rate,
indicative of the nuclear relaxation due to scattering of magnons [30].
According to Beeman and Pincus [31], in the AF state for magnetic insulators,
1/$T_{1}$ is mainly driven by magnon processes, leading to a strong
temperature dependence of 1/$T_{1}$. Here, we calculate the temperature
dependence of 1/$T_{1}$ based on the two-magnon (Raman) process. With a long
wave approximation for magnon dispersion relation, 1/$T_{1}$ can be given by
[31]

$\frac{1}{T_{1}}=CT^{3}\int_{T_{AE}/T}^{\infty}\frac{x}{e^{x}-1}dx.$ (10)

Here, the fitting parameters are the proportional constant $C$ and $T_{\rm
AE}$. $T_{\rm AE}$ is the anisotropy gap in the spin wave spectrum. The solid
red curve in Fig. 7 shows the calculated results with $C$ = 2.5 (1/sK3) and
$T_{\rm AE}$ = 5.7 K, which reproduced the experimental data well.

We also calculated the temperature dependence of 1/$T_{1}$ for the two magnon
process discussed in triangular antiferromagnets CsNiBr3 [32] and NiGa2S4
[33], and a Kagome lattice antiferromagnet KFe3(OH)6(SO4)2 [34]. In this case,
the temperature dependence of 1/$T_{1}$ is given by

$\frac{1}{T_{1}}=C^{\prime}{}T^{5}\int_{T_{\rm AE}/T}^{T_{\rm
m}/T}\left[x^{2}-\left(\frac{T_{\rm
AE}}{T}\right)^{2}\right]\left[x^{2}+\left(\frac{T_{m}}{T}\right)^{2}\right]\frac{e^{x}}{(e^{x}-1)^{2}}dx,$
(11)

where $T_{m}$ is defined by $\hbar\omega_{\rm m}/k_{\rm B}$ with the maximum
frequency of the spin wave ($\omega_{\rm m}$), and $T_{\rm m}=zSJ$ ($\sim$ 29
K) where $z=6$ is the number of nearest neighbor of Fe ions surrounding a Fe
ion. As shown by the green dotted curve, which is the calculated result with
$C^{\prime}$ = 1.3$\times$10-3 (1/sK5) and $T_{\rm AE}$ = 5.7 K, the model
also reasonably reproduces experimental results. Since both models reproduce
the experimental data very well, the magnetic fluctuations in the AFM state of
NaFe3(PO4)2 are considered to be well described by the conventional magnon
scattering without showing any additional magnetic fluctuations due to spin
frustrations, if any.

The magnon picture also well reproduces the temperature dependence of the
sublattice magnetization at low temperatures. Using the long wave
approximation, the temperature dependence of the sublattice magnetization is
given by

$M(0)-M(T)\propto T^{2}\int_{T_{AE}/T}^{\infty}\frac{x}{e^{x}-1}dx.$ (12)

The red curve in Fig. 5 shows the calculated results with the same value of
the parameter ($T_{\rm AE}$ = 5.7 K), which can reproduce the temperature
dependence of $H_{\rm int}$. The result again shows the conventional nature of
the AFM state in NaFe3(PO$4$)2.

Finally, we briefly mention magnetic anisotropy. Experimentally the magnetic
anisotropy in NaFe3(PO$4$)2 is clearly observed in the present $T_{1}$ data in
the AFM state and also the magnetization measurements reported previously
[11]. Since the orbital moments of Fe3+ ions are zero in the ${}^{6}S$ ground
state for the 3$d^{5}$ electron configuration, one may expect no single-ion
anisotropy in the first order. However, because of the second-order
perturbation between the ground and excited states, the single ion type
anisotropy can be induced by spin-orbit interactions. Dzyaloshinsky-Moriya
interaction also may produce magnetic anisotropy in the magnetically ordered
state.

###  3.5 Nuclear spin-spin relaxation rate ($1/T_{2}$)

Figure 9: (a) 31P spin-spin relaxation rate $1/T_{2}$ vs. $T$ measured at 32.5
MHz. (b) Temperature dependence of $\alpha^{\prime}$ [$\equiv$
(1/$T_{1}$)/(1/$T_{2}$)]. (c) 1/$T_{1}$ vs. $\varepsilon$ and
$\alpha^{\prime}{}$/$T_{2}$ vs $\varepsilon$ where $\varepsilon$ = $(T-T_{\rm
N})/T_{\rm N}$ with $T_{\rm N}$ = 11 K.

The temperature dependence of $1/T_{2}$ measured at $f$ = 32.5 MHz is shown in
Fig. 9(a), and the trend agrees with that of $1/T_{1}$ behavior. It can be
seen that for $T\leq 100$ K, $1/T_{2}$ gradually decreases down to $T\simeq
15$ K. Then 1/$T_{2}$ increases sharply and exhibits a peak at $T_{\rm N}\sim
11$ K, a clear indication of AFM long-range order. Below $T_{\rm N}$,
$1/T_{2}$ decreases and levels off at low temperatures. It is noted that the
values of 1/$T_{2}$ are slightly greater than those of $1/T_{1}$ measured at
the same magnetic field above $T_{\rm N}$, although there is a large
difference between them in the AFM ordered state.

In general, 1/$T_{2}$ is related to 1/$T_{1}$ and can be written for the case
of $I$ = 1/2 as

$\frac{1}{T_{2}}=\left(\frac{1}{T_{2}}\right)^{*}+\left(\frac{1}{2T_{1}}\right)+F_{z}(0),$
(13)

with

$\frac{1}{T_{1}}=F_{\perp}(\omega_{\rm N}).$ (14)

Here, $F_{\alpha}(\omega)$ is the spectral density of the longitudinal
($\alpha$ = $z$) and transverse ($\alpha$ = $\perp$) components of the
fluctuating local field ($h_{\alpha}$), and is described by

$\displaystyle F_{\alpha}(\omega_{\rm N})=\frac{1}{2}\gamma_{\rm
N}^{2}\int_{-\infty}^{+\infty}\langle\ h_{\alpha}(t)h_{\alpha}(0)\ \rangle{\rm
exp}(i\omega_{\rm N}t)dt.$ (15)

$(1/T_{2})^{*}$ is due to the nuclear dipole-dipole interaction [36, 35] and
is independent of temperature. This can be attributed to the nearly
temperature-independent 1/$T_{2}$ observed at low temperatures in the AFM
state. The lowest value of 1/$T_{2}$ $\sim$ 1300 (1/s) corresponds to the
spectral width $0.8$ Oe, which will be of the order of nuclear dipolar field
at the P site.

Thus the temperature dependent part of 1/$T_{2}$ comes from the other two
terms originating from $F_{z}(0)$ and $F_{\perp}(\omega_{\rm N})$. To see how
those contributions to 1/$T_{2}$ change with temperature, we plotted
$\alpha^{\prime}$ [$\equiv$ (1/$T_{1}$)/(1/$T_{2}$)] shown in Fig. 9(b) where
we used the 1/$T_{1}$ data measured at the same frequency of 32.5 MHz.
$\alpha^{\prime}$ is nearly constant above 50 K, and starts increasing below
50 K with decreasing temperature and becomes close to $\alpha^{\prime}$ $\sim$
0.9 near $T_{\rm N}$. Just below $T_{\rm N}$, $\alpha^{\prime}$ shows a sudden
decrease due to AFM ordering. Using $\alpha^{\prime}$ = 0.6 above 50 K,
$F_{z}(0)$ is estimated to be 1.17$F_{\perp}(\omega_{\rm N})$. This suggests
that $F_{z}(0)$ and $F_{\perp}(\omega_{\rm N})$ contribute nearly equally to
1/$T_{2}$. Below 50 K, on the other hand, the increase in $\alpha^{\prime}$
suggests that the contribution of $F_{z}(0)$ to 1/$T_{2}$ becomes small. If we
take $\alpha^{\prime}$ = 0.9, $F_{z}(0)$ is calculated to be
0.61$F_{\perp}(\omega_{\rm N})$. As we observed the development of AFM spin
fluctuations below 50 K, the increase of $\alpha^{\prime}$ could be due to the
AFM spin fluctuations. Since $F_{\perp}(\omega_{\rm N})$, that is, 1/$T_{1}$,
picks up the hyperfine fluctuations at NMR frequency of $\omega_{\rm N}$ while
$F_{z}(0)$ relates to the hyperfine fluctuations at $\omega$ = 0, the increase
of $\alpha^{\prime}$ may suggest that the spectrum density of the AFM spin
fluctuation does not extend to $\omega$ = 0 in the paramagnetic state.

On the other hand, we observed the clear divergence behaviors in 1/$T_{2}$ and
$1/T_{1}$ very near $T_{\rm N}$ which is due to the critical slowing down of
magnetic fluctuations to $\omega$ = 0, characteristic of second order phase
transition. To see the critical behavior we plotted 1/$T_{1}$ and
$\alpha^{\prime}$/$T_{2}$ near $T_{\rm N}$ in Fig. 9(c) as a function of
$\varepsilon$ = $(T-T_{\rm N})/T_{\rm N}$. The divergent behaviors are fitted
with the power law of $1/T_{1}$ (or $\alpha$/$T_{2}$) $\propto$
$\varepsilon^{-0.38}$. This value of power is close to 0.3 expected for the 3D
Heisenberg spin system but far from 0.8 for the 2D Heisenberg spin system
[37], indicating that the AFM transition is most likely characterized by the
3D nature.

##  4 Conclusion

The dynamic and static properties of the 2D anisotropic triangular lattice
compound Na3Fe(PO${}_{4})_{2}$ have been investigated from a microscopic point
of view by the NMR technique. The 31P NMR spectra, nuclear spin-lattice
relaxation rate 1/$T_{1}$ and nuclear spin-spin relaxation rate 1/$T_{2}$
results show the antiferromagnetic ground state below $T_{\rm N}\sim$ 11 K. In
the paramagnetic state above $T_{\rm N}$, from the analysis of NMR shift
$K(T)$ using the high-temperature series expansion, $J/k_{\rm B}$ was
estimated to be $\simeq 1.9$ K, which is consistent with the previously
reported value of $J/k_{\rm B}\simeq 1.8$ K. The antiferromagnetic spin
fluctuations were found to be developed below $\sim$50 K from the temperature
dependence of $1/T_{1}T\chi$. In the magnetically ordered state below $T_{\rm
N}$, from the observation of the characteristic rectangular shape of the NMR
spectra, the AFM state was found to be a commensurate AFM state in its ground
state, consistent with the neutron diffraction measurements [11]. The static
and dynamic magnetic properties in the AFM state were revealed to be well
described by the conventional spin-wave theory for a 3D antiferromagnet with a
spin-anisotropy energy gap of 5.7 K. Our NMR results indicate that the
magnetically ordered state of Na3Fe(PO4)2 is a conventional 3D
antiferromagnetic state without any obvious effects of spin frustrations in
its ground state.

##  5 Acknowledgments

The research was supported by the U.S. Department of Energy (DOE), Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering. Ames
Laboratory is operated for the U.S. DOE by Iowa State University under
Contract No. DE-AC02-07CH11358. SJS and RN would like to acknowledge SERB,
India, for financial support bearing sanction Grant No. CRG/2019/000960. SJS
is supported by the Prime Minister’s Research Fellowship (PMRF) scheme,
Government of India.

##  6 References

## References

  * [1] Li Y, Gegenwart P and Tsirlin A. A 2020 J. Phys.: Condens. Matter 32 224004
  * [2] Anderson P 1973 Mater. Res. Bull 8 153
  * [3] Huse D A and Elser V 1988 Phys. Rev. Lett 60 2531\.
  * [4] Jolicoeur T and Le Guillou J C 1989 Phys. Rev. B 40 2727
  * [5] Sen A, Wang F, Damle K and Moessner R 2009 Phys. Rev. Lett. 102, 227001
  * [6] Tapp, J, dela Cruz C R, Bratsch M Amuneke, N E, Postulka L, Wolf B, Lang M, Jeschke H O, Valentí R, Lemmens P and Möller A 2017 Phys. Rev. B 96 064404
  * [7] Sebastian S J, Somesh K, Nandi M, Ahmed N, Bag P, Baenitz M, Koo B, Sichelschmidt J, Tsirlin A A, Furukawa Y and Nath R 2021 Phys. Rev. B 103 064413
  * [8] Amuneke N E, Tapp J, de la Cruz C R and Möller A 2014 Chem. Mater. 26 5930
  * [9] Belkhiria M, Laaribi S, Amara A and Amara M B1998 Ann. Chim.: Sci. Mater. 23 117
  * [10] Morozov V, Lazoryak B, Malakho A, Pokholok K, Polyakov S and Terekhina T 2001 J. Solid State Chem. 160 377
  * [11] Sebastian S J, Islam S S, Jain A, Yusuf S M, Uhlarz M and Nath R 2022 Phys. Rev. B 105 104425
  * [12] Carter G S, Bennett L H and Kahan D J, 1977 Metallic Shifts in NMR (Pergamon Press)
  * [13] Delmas C, Le Flem G, Fouassier C and Hagenmuller P 1978 J. Phys. Chem. Solids 39 55
  * [14] Schmidt H J, Lohmann A and Richter J, 2011 Phys. Rev. B 84 104443
  * [15] Yamada Y and Sakata A 1986 J. Phys. Soc. Jpn. 55 1751
  * [16] Ranjith K. M, Majumder M, Baenitz M, Tsirlin A A and Nath R 2015 Phys. Rev. B 92 024422
  * [17] Ranjith K M, Nath R, Majumder M, Kasinathan D, Skoulatos M, Keller L, Skourski Y, Baenitz M and Tsirlin A A 2016 Phys. Rev. B 94 014415
  * [18] Nath R, Ranjith K M, Roy B, Johnston D C, Furukawa Y and Tsirlin A A 2014 Phys. Rev. B 90 024431
  * [19] Bramwell S T and Holdsworth P C W 1994 Phys. Rev. B 49 8811
  * [20] Collins M F 1989 Magnetic Critical Scattering (Oxford University Press, New York)
  * [21] Ozeki Y and Ito N, 2007 J. Phys. A: Math. Theor. 40 R149
  * [22] Nath R, Furukawa Y, Borsa F, Kaul, E E, Baenitz M, Geibel C, Johnston D C, 2009 Phys. Rev. B 80 214430
  * [23] Sebastian S J, Somesh K, Nandi M, Ahmed N, Bag P, Baenitz M, Koo B, Sichelschmidt J, Tsirlin A A, Furukawa Y and Nath R 2021 Phys. Rev. B 103, 064413
  * [24] Takigawa M, Motoyama N, Eisaki H and Uchida S 1996 Phys. Rev. Lett. 76 4612
  * [25] Hone D, Scherer C and Borsa F 1974 Phys. Rev. B 9 965
  * [26] Ajiro Y, Nakajima Y, Furukawa Y and Kiriyama H 1978 J. Phys. Soc. Jpn. 44 420
  * [27] Furukawa Y, Iwai A, Kumagai K and Yakubovsky A 1996 J. Phys. Soc. Jpn. 65 2393
  * [28] Yogi A, Ahmed N, Nath R, Tsirlin A A, Kundu S, Mahajan A V, Sichelschmidt J, Roy B, and Furukawa Y 2015 Phys. Rev. B 91, 024413 (2015).
  * [29] Smerald A and Shannon N 2011 Phys. Rev. B 84 184437
  * [30] Belesi M, Borsa F and Powell A K 2006 Phys. Rev. B74 184408
  * [31] Beeman D and Pincus P 1968 Phys. Rev. 166 359
  * [32] Maegawa S 1995 Phys. Rev. B 51 15979
  * [33] Takeya H, Ishida K, Kitagawa K, Ihara Y, Onuma K, Maeno Y, Nambu Y, Nakatsuji D, MacLaughlin D E, Koda A and Kadono R 2008 Phys. Rev. B 77 054429
  * [34] Nishiyama M, Maegawa S, Inam T and Oka Y 2003 Phys. Rev. B 67 224435
  * [35] Slichter C P 1989 $Principles$ $of$ $Magnetic$ $Resonance$ (Berlin: Springer)
  * [36] Abragam A 1961 $The~{}Principles~{}of~{}Nuclear~{}Magnetism$ (Oxford: Clarendon Press)
  * [37] Benner H and Boucher J P 1990 Magnetic Properties of Layered Transition Metals Compounds (Dordrecht: Kluwer Academic Publishers)


# PE-GAN: Prior Embedding GAN for PXD Images at Belle II

Hosein Hashemi
Faculty of Physics
Ludwig Maximilians University of Munich,
Munich, Germany
<EMAIL_ADDRESS>
Nikolai Hartmann
Faculty of Physics
Ludwig Maximilians University of Munich,
Munich, Germany
<EMAIL_ADDRESS>
Thomas Kuhr
Faculty of Physics
Ludwig Maximilians University of Munich,
Munich, Germany
<EMAIL_ADDRESS>
Martin Ritter
Faculty of Physics
Ludwig Maximilians University of Munich,
Munich, Germany
<EMAIL_ADDRESS>
Matej Srebre
Faculty of Physics
Ludwig Maximilians University of Munich,
Munich, Germany
<EMAIL_ADDRESS>

###### Abstract

The pixel vertex detector (PXD) is an essential part of the Belle II detector
recording particle positions. Data from the PXD and other sensors allow us to
reconstruct particle tracks and decay vertices. The effect of background hits
on track reconstruction is simulated by adding measured or simulated
background hit patterns to the hits produced by simulated signal particles.
This model requires a large set of statistically independent PXD background
noise samples to avoid a systematic bias of reconstructed tracks. However,
data from the fine-grained PXD requires a substantial amount of storage. As an
efficient way of producing background noise, we explore the idea of an on-
demand PXD background generator using conditional Generative Adversarial
Networks (GANs), adapted by the number of PXD sensors in order to both
increase the image fidelity and produce sensor-dependent PXD hitmaps.

_Keywords_ GANs $\cdot$ High-granularity detector simulation $\cdot$ PXD Belle
II $\cdot$ Contrastive Learning

##  1 Introduction

Belle II [1] is a B-factory experiment located at the SuperKEKB [2]
accelerator in Tsukuba, Japan operating at peak luminosity of $2\times
10^{34}$ cm-2s-1. The accelerator collides electrons and positrons at a center
of mass energy of $\sqrt{s}\approx 10.6$ GeV, which prompts to $\Upsilon(4S)$
resonance state that decays almost exclusively to entangled pairs of either
charged or neutral $B$ mesons. These processes are called signal decays which
are isotropic events. We also have other processes which are jet-like
continuum events $(e^{+}e^{-}\longrightarrow q\bar{q},~{}q=u,d,s,c)$. The
precise measurement of $CP$-violation and the search for rare or “forbidden”
decays of the B-meson and the tau-lepton as signals for New Physics relies
heavily on a large number of recorded events and the precision with which
B-meson and lepton decay vertices can be reconstructed. In order to do so,
Belle II plans to collect a 50 times larger dataset than its predecessor,
Belle. The accelerator’s upgrade resulted in the desired increase of the
collision rate of particles. Thus, we expect larger background at design
luminosity for all sub-detectors. To profit from the much larger recorded
dataset, the systematic uncertainties of simulated datasets must be well under
control. This requires an accurate incorporation of beam background effects in
the detector simulation.

Various background processes lead to signals in the detector that do not
originate from the physics processes of interest. These processes can be
categorised into two types, beam-induced and luminosity dependent processes.
The beam-induced processes originate from the synchroton radiation and
collisions of beam particles with residual gas in the beampipe, bending
magnets or particles within a bunch, whereas luminosity dependent processes
comprise electron-positron collisions leading to “non-interesting” physics
processes such as Bhabha scattering or two-photon processes.

There are two methods to include and simulate these beam background effects,
the background mixing method and the background overlay method. In the
background mixing method the pre-simulated Geant4 [3] background hits are
added to simulated energy deposit as SimHits from a signal event. Then, both
background and signal contributions are merged and digitized at the same time.
This has the advantage that the subsequent simulation of the digitization
correctly takes into account cases where energy deposits from background and
signal particles contribute to the same channel. On the other hand the
comparison of background simulations using this so-called mixing technique
with measurements show sizable discrepancies. The background overlay method
uses random trigger events along with real data to create a background event
sample. This works on digit level and thus needs a detector-specific treatment
of cases where background and signal contribute to the same channel.
Nevertheless it is expected to yield a more realistic simulation of beam
background effects.

##  2 Pixel Vertex Detector

The Pixel Vertex Detector (PXD) [4] was designed to be the innermost semi-
conductor sub-detector of Belle II. The PXD measures the position of
traversing particles originating from a particle collision in order to do
precise reconstruction of decay vertices. As a result it is placed very
closely to the interaction point, hence exposed to a huge amount of background
radiation. The designed configuration of PXD consists of 40 modules within 2
layers of detector. The inner layer has 16 modules assembled in 8 planar
ladders. The outer layer is composed from 12 ladders or 24 modules. Each
module consists of a $250\times 768$ pixel matrix of DEPFET [5] pixel sensors
which culminates in a $250\times 768$ hitmap. Currently only 20 modules are
installed. With a readout time of $20~{}\mu s$ the PXD integrates over about
$5000$ bunch crossings. Thus, the majority of pixels with a charge above
threshold stems from the beam background. Pixel signals are only read out if
their value exceeds a threshold. In our simulated data the threshold is at 7
ADU. We will address this in section 4.2.

###  2.1 The Problem

A problem with the background overlay method is that large amounts of
resources are required for the readout, storage, and distribution of the
random trigger events. A single event in the format to be used for the
background overlay has a size of about $200$ kB. This is roughly a factor 20
more than the size of events in the format used for analyses. Furthermore, the
beam background files have to be distributed to all sites where run-dependent
Monte-Carlo samples are produced. This multiplies the storage demands by the
number of sites. In particular, smaller sites do not have the storage
resources with the required size and bandwidth to host full-beam background
datasets. Even if one copies the background files from the main server, and
wants to access them through the network, each job has to access multiple
files, for example for producing $1$ Tb of MC files, one has to copy $20$ Tb
PXD background files which dramatically increases the data transfer and
networking traffic.

Our solution to this problem is to make it an end-to-end procedure. First, to
generate these hitmaps on the fly and upon request for each sensor with a
conditional generative agent, then extract PXD hits data from the fired
pixels. The only storage cost would be the weights of the generative model.

##  3 PXD Background Generation

This work is in the direction of the "Fast Simulation" campaign in Particle
Physics that had been started after the introduction of the Generative
Adversarial Networks (GAN), by Goodfellow et al. [6] for image data. GAN is a
flexible unsupervised deep learning architecture that can be trained to learn
a function to match a given data distribution. GANs have been widely used in
Particle Physics in order to have a faster detector simulation for LHC [7]
[8], but it has not been used for Pixel detector. GANs involve two networks,
the Generator, and the Discriminator, in a two-player zero-sum game whose
purpose, respectively, is to map random noise to samples and discriminate real
and generated samples. In the end, the goal is to find a Nash equilibrium to a
two-player min-max problem. In our case, we need a model to generate sensor-
dependent PXD background hitmaps (images) as fast as possible with the highest
fidelity. For the moment with the current approach, we ignore correlations
between the two layers, since most of the background is coming from out-of-
time events. As we usually do not have the corresponding SVD hits for two
correlated PXD background hits caused by the same particles, so no track would
be reconstructed. Thus, the expected effect on tracking performance from
neglecting the correlations is small but still has to be studied in detail.

Our data in each event consists of 40 sensor-dependent grey-scale $250\times
768$ images. For training, we used 40000 MC simulated events. And for
evaluation, we have 10000 MC simulated events. The events are simulated for
two times the expected background level at design luminosity. We are also
planning to use the real data (random trigger) in the future as well.

Here we are facing several challenges:

  1. 1.

Having high-resolution images. In general with GANs, there is no guarantee of
balance between the training of the generator and discriminator. As a
consequence, one network may be more powerful than the other, which in most
cases is the discriminator. When the discriminator becomes much stronger than
the generator, the fake images become too easy to be separated from real ones,
thus reaching a point where gradients from the discriminator vanish, providing
no guidance for further training of the generator. This happens more
frequently when one is dealing with high-resolution images due to the
difficulty of generating meaningful high-frequency details. On the other hand,
having less powerful discriminator results in a mode collapse, as with the
generator greedily optimizing its loss function in order to “fool” the
discriminator, each generator update is a partial collapse towards mode
collapse [9]. A crucial aspect in maintaining a balance between solutions for
vanishing Gradients and mode collapse is a discriminator with high
generalization capabilities. Hence, in the end, it is provided by deeper
networks with non-local blocks along with having PXD ladder numbers (sensor
IDs) as conditional classes in an attempt to capture global and local
information and dependencies which are notoriously memory consuming.

  2. 2.

Having sparse images. Convolutional Neural Network (CNN) as a part of the
discriminator has trouble with sparse data [10]. They require dense data to
learn well. Naively applying CNNs to sparse data only works partially as they
are sensitive to missing data. Thus, we need more sophisticated methods and
prior information such as sensor IDs to be injected into the learning
procedure.

  3. 3.

The lack of evaluation while training. When one deals with the generic image
dataset such as ImageNet [11], one can use the Inception score (IS) [13] or
the Frechet Inception Score (FID) [14] for training evaluation and validation
to choose the best model as the GAN training is utterly unstable. But in our
case, due to the specificity of our dataset, the commonly used evaluation
metrics for natural images don’t apply.

###  3.1 The Model

After considering a mixture of the models and structures, and tuning many
hyperparameters by trial and error minimization based on our evaluation
methods mentioned in section 5, we chose a modified version of BigGAN-deep
[12] with contrastive conditioning [16] and consistency regularization [18].

We started with adding class labels which were 40 sensor numbers to the model
in an extra embedding layer to our former model, WGAN-gp [20], which acted as
a baseline to our study. The problem was slow event-based image generation in
comparison to the overall simulation time due to its large kernel filter size,
and not capturing entirely the pixel intensity distribution. Then, we went
through more complex but faster models due to the lower kernel size and
special residual block which were designed for high-resolution image
generation such as BigGAN [12]. During the training, our models undergo
training collapse, thus early stopping was crucial. As suggested in BigGAN’s
paper [12], we monitored and shown in Fig. 2, the first singular values of
each weight matrix for both the generator and the discriminator grow
throughout training, but whenever we had mode collapse or discriminator
overfitting to the training set, the first singular values explodes
drastically for the generator. We tried many treatments, such as lowering the
learning rate for the discriminator, using different regularization methods
such as Orthogonal Regularization, Minibatch Discrimination [13], Dropout, and
L2 regularization, Latent Optimization [21], Feature Quantisation [22] did not
work either. As a result, we thought that we have to go for a simpler model.

Then, based on the context of our image dataset, and memory limits, along with
the feedbacks and failures from the model during the training, we decided to
use a simpler model. Thus, we used SAGAN [15]. Self-Attention GAN is an
architecture that allows both the generator and discriminator to model long-
range dependencies over the image manifold. The key idea is to make the
generator able to "see" not only the convolutional kernel range but also the
global details. Normally, Convolutions as local blocks have very nice
properties such as parameter sharing and translation invariance and they learn
in a hierarchical order. Nevertheless, when it comes to complex structures and
elaborated geometrical contours, it requires long-range details. This is where
attention enters. With this choice of architecture, we learned a lot about how
the difference between the models would affect our result and the final
evaluation result was surpassed by a deeper model with fewer intricacies with
respect to BigGAN but more layers than SAGAN.

In our model, we increased the depth of each block in both the generator and
the discriminator, with 32 multiplication channels in each layer. After some
training and ablation study over BigGAN-deep, we decided to remove the
Projection Discrimination which was originally introduced in SNGAN [23].
Projection Discrimination’s objective is to maximize the projection of the
embeddings of real images on the corresponding target embeddings while
minimizing the inner-product values when the images are fake. However, since
it only utilizes data-to-class relation, it results in the discriminator’s
overfitting and training collapse problems [12]. Hence, we replaced it with
the contrastive conditioning discrimination mechanism [16]. It is based on a
conditional contrastive loss that comes from metric learning that considers
data-to-data relations in each batch. It pulls the multiple-image embeddings
closer to each other in the embedding space when the class labels are the
same. In other words, the discriminator can capture not only data-to-class but
also data-to-data relations between samples.

Within the context of metric learning loss, due to the fact that we are
dealing with a close-conditional image generation task, we also added an
embedding of the prior information along with the class labels to both the
discriminator and generator. The relational embedding module in the
discriminator consists of a linear projection of semantic information along
with a self-attention over the concatenated version of label embeddings and
semantic information projection. Let $\mathbf{X}=\\{x_{1},...,x_{m}\\}$, where
$x_{i}\in\mathbb{R}^{W\times H}$ be sampled training images and
$\mathbf{y}=\\{y_{1},...,y_{m}\\}$, where $y_{i}\in\mathbf{Z}$ be the
corresponding labels from 1 to 40, and $\mathbf{s}=\\{s_{1},...,s_{m}\\}$,
where $s_{i}\in\mathbf{R}$ be the mean occupancies for each class of sensors.
Also, we define a linear projection function
$l^{(d)}:\mathbb{R}^{k}\rightarrow\mathbb{S}^{d}$ and an embedding function
$\mathbf{e}:y_{i}\rightarrow\mathbb{R}^{d/4}$. For the relational embedding,
we have 40 number of tokens of dimension $d$ as query, key, and value from the
output of the concatenation, and a dot-product self-attention [17] which
calculates the weighted average of feature representations with the weight
proportional to the dot-product similarity score between pairs of
representation,

$Q=l^{(d)}~{}[l^{(d/4)}(\mathbf{s})~{}||~{}\mathbf{e}]$
$A=softmax~{}(\frac{QK^{T}}{\sqrt{d}})~{}V.$

Then, using the linear projection of the output of the discriminator, one can
define the conditional contrastive loss are as,

$L(x_{i},y_{i},s_{i};t)=-\log\left(\frac{\exp{(h_{i}^{T}a_{i}/t)}+\sum^{m}_{k=1}\mathbf{1}_{k=i}.\exp{(h_{i}^{T}h_{k}/t)}}{\exp{(h_{i}^{T}a_{i}/t)}+\sum^{m}_{k=1}\mathbf{1}_{k\neq
i}.\exp{(h_{i}^{T}h_{k}/t)}}\right)$ $H=l^{(d)}[D(x_{i})]$

Where $a_{i}\in A$ and $h_{i}\in H$ . Minimizing the above contrastive loss
will reduce distances between the embeddings of images with the same labels
while maximizing the others. It considers the data-to-data relations
$h_{i}^{T}h_{k}$ and data-to-class relations $h_{i}^{T}a_{i}$ without
augmentations. Thus, we help the discriminator to be guided by not only the
image but also by the fine-grained prior information as an inductive bias
injection. This prior semantic information is based on different mean
occupancies from the simulated images. This extra information injection
boosted both the performance regarding the mean-occupancy distribution and the
stability of the training for almost 20000 more iterations.

We also used consistency regularization [18] within our training which
prolonged the training stability significantly, and delayed collapse. The idea
behind consistency regularization is to penalize the sensitivity of the
classifier to some semantics-preserving noise. In our case, the noise
perturbation is image flipping and random shift as images with these augmented
artifacts exist in the generator’s outputs.

Eventually, our model, as depicted in Table 1, using the above technologies
along with orthogonal initialized [24] residual blocks [12] for both the
generator and the discriminator and Hinge loss [15], succeeded in the sensor-
dependent PXD hitmap generation. We trained our model with the PyTorch
framework [19] over Tesla V100 GPU. This model generates 40 hitmaps for an
event in around 6 seconds on one CPU.

$z\in\mathbb{R}^{128}\sim\mathcal{N}(0,I)$
---
Embed(y) $\in\mathbb{R}^{128}$
$ch=32$
Linear(128+128) $\rightarrow 4\times 12\times 16ch$
ResBlock $16ch\rightarrow 16ch$
ResBlock up $16ch\rightarrow 16ch$
ResBlock $16ch\rightarrow 16ch$
ResBlock up $16ch\rightarrow 8ch$
ResBlock $8ch\rightarrow 8ch$
ResBlock up $8ch\rightarrow 8ch$
ResBlock $8ch\rightarrow 8ch$
ResBlock up $8ch\rightarrow 4ch$
Non-Local Attention Block
ResBlock $4ch\rightarrow 4ch$
ResBlock up $4ch\rightarrow 2ch$
ResBlock $2ch\rightarrow 2ch$
ResBlock up $2ch\rightarrow ch$
BN, ReLU, $3\times 3$ Conv $ch\rightarrow 1$
Tanh
(a) Generator

zero padded gray-scale image
---
$x\in\mathbb{R}^{256\times 768\times 1}$
ConvBlock $1\rightarrow ch$
ResBlock down $ch\rightarrow 2ch$
ResBlock $2ch\rightarrow 2ch$
ResBlock down $2ch\rightarrow 4ch$
ResBlock $4ch\rightarrow 4ch$
ResBlock down $4ch\rightarrow 8ch$
ResBlock $8ch\rightarrow 8ch$
ResBlock down $8ch\rightarrow 8ch$
ResBlock $8ch\rightarrow 8ch$
ResBlock down $8ch\rightarrow 16ch$
ResBlock $16ch\rightarrow 16ch$
ResBlock down $16ch\rightarrow 16ch$
ResBlock $16ch\rightarrow 16ch$
ReLU, Global sum pooling
Self-Attention [Proj(Embed(y) || Linear(s))],
Norm(Linear512$\rightarrow 1024$)
Conditional Contrastive loss
(b) Discriminator

Table 1: The Architecture of the Generator and the Discriminator. "$ch$"
denotes the number of multiplicative channels, "ResBlock*", the residual
blocks introduced in BigGAN-deep, "BN", the batch normalization, "ReLU", the
rectified linear unit, "Norm", a normalization operator, "Proj" and "Linear",
a linear projection or an MLP.

##  4 Model Validation

For model evaluation we have two categories of methods. The ones that we use
within the training process, and the ones after training for validation.

  1. 1.

During training:

     * •

Image visualization (Turing Test)

     * •

Monitoring singular values for Weight matrices

  2. 2.

After training:

     * •

Image pixel intensity analysis

     * •

Sensor number occupancy distribution

     * •

Physics analysis over impact parameters resolutions

###  4.1 Training Evaluation Methods

Fig.1 shows a comparison between a generated hitmap and a real one for a
sensor in the inner layer. One can generate sensor-dependent images for any
layer and ladder number. It is very important for the model to pass the first
qualitative evaluation test. Blurriness, image fidelity, general clustering
details, and sharpness are some visual factors that one can check. Needless to
mention that, we are also developing an automated and quantitative version of
this evaluation, discussed in the Outlook section.

(a) a

(b) b

Figure 1: Comparison of color-reversed PXD image (up) and generated image
(down)

Regarding the first singular value of the weight matrix for different layers,
our observation was in accord with [12] that when the loss for the
discriminator approaches zero constantly during training, it is due to the
collapse. However, we also observed that for the discriminator, the spectra
are very noisy, and when either mode collapse or discriminator memorization
happens, the singular values drop suddenly. The most important thing to
observe is the absolute change. In Fig.2, one can see the first singular
values of the weight matrices of each layer of the model. They can be
efficiently computed using the Arnoldi iteration method [25]. The first large
drop and consequently collapse occurs around iteration 25000 with the green
and orange layer which can be seen in Fig.2a. As a result, the image contents
gradually become complete nonsense. In order to address this issue and delay
the collapse we tried and evaluated many methods such as Dropout, Orthogonal
Regularization for both the generator and discriminator, Spectral
regularization [26] and also gradient penalty [27]. However, the gained
performance was not desirable in terms of stability endurance. Then, after
modifying the BigGAN-deep model with contrastive conditioning over the
discriminator, we could achieve proper stability in the end.

(a) Collapsed model (b) Our Model

Figure 2: First singular values of the discriminator for each layer in (a) a
model with collapse (b) our model. The most important thing to notice here is
the sudden drop in the values over a absolute large range which is an
indication of discriminator’s overfitting.

###  4.2 Validation Methods

For the post-training validation process, we used 10000 new simulated events.
First, we observed the pixel intensity values. Pixel value $0$ means complete
blackness. Based on the MC simulation, there is a gap between 0 and 7 ADU as
in fig.3. Below pixel intensity $7$, we used a function to cut off the
residual pixel values. As a result, we will recover the separation between the
two parts of the diagram. However, we will lose many pixel values which in
principle correspond to a specific background process in the image. Moreover,
the non-linearity of the intensity curve is not fully learned yet. We will
address this issue in the discussion section.

Figure 3: Pixel intensity value comparison. Obviously a large portion of the
images are totally dark. There is a very narrow discrepancy between the Geant4
simulated and GAN generated pixel values.

Our other metric for validation is the mean occupancy frequency comparison
between the MC simulated test set and the generated images in Fig.4. This is a
consequence of considering sensor IDs in our class-conditional GAN model. As a
result of considering the labels, the model partially learned the frequency of
PXD hits over each sensor which is of real importance. Subsequently, this
indicates the generation of hitmaps closer to reality with better diversity
and fidelity. The discrepancy between the real and generated mean occupancy
over the hitmaps can be lowered by more tuning over the model. We propose more
solutions in sec. 5.

Figure 4: Sensor dependent mean occupancy comparison between Geant4 simulated
and GAN generated images

In the end, our last metric is a quantifiable comparison between the generated
hitmaps and the simulated ones by studying the impact parameter resolutions.
To do so, we simulated $\Upsilon(4S)$ signal events and considered three
background overlay scenarios to compare with each other:

  * •

No background in the PXD

  * •

Geant4 simulated PXD background

  * •

Generated PXD background with modified BigGAN-deep

Figure 5: Comparison of the resolution of the signed distance in the
transverse plane from the pivotal point to the helix, between simulated and
generated background and no background

The main effect of the background in the PXD is that wrong hits are attached
to tracks. As a result, the effect of the background is visible over the tails
of the distribution. The difference between the Geant4 simulated and
conditional GAN is marginal and our model successfully reproduce Geant4
simulated background effects. In order to quantify the difference between
resolutions, we introduce the 2-point Wasserstein distance (Frechet Distance)
[28] which is a similarity measure for Gaussian distributions. For two
multivariate Gaussian distributions with means $\mu_{x}$ and $\mu_{y}$ and
covariance matrices, $\Sigma_{x}$ and $\Sigma_{y}$, the Frechet distance is
defined as

$W^{2}(x,y)=|\mu_{x}-\mu_{y}|^{2}+\operatorname{Tr}(\Sigma_{x}+\Sigma_{y}-2\sqrt{(\Sigma_{x}\Sigma_{y})}~{}).$
(1)

Looking at the two-point Wasserstein distance function between the
distributions, the lower the Frechet Distance (FD), the better the image
quality from the physics point of view. As an outcome of conditioning over the
sensor IDs, the Frechet Distance score became even smaller in comparison to
the former model.

Model | Geant4. vs None. | GAN. vs None. | GAN. vs Geant4.
---|---|---|---
modified BigGAN-deep | $5.46\mathrm{e}{-3}$ | $\mathbf{6.70\mathbf{e}{-3}}$ | $\mathbf{2.51\mathbf{e}{-3}}$
WGAN-gp | $5.46\mathrm{e}{-3}$ | $6.98\mathrm{e}{-3}$ | $2.82\mathrm{e}{-3}$
Table 2: FD scores comparison between models in three background scenarios.

##  5 Summary and outlook

We developed a class-conditional GAN, based on PXD images and their sensor ID.
As a result, we could achieve more reasonable agreement with real data and now
are able to generate sensor-dependent PXD background hitmaps. Using GAN to
generate the background data, we will reduce the storage consumption with the
cost of additional CPU resources. In our path, we faced many challenges such
as severe mode-collapse and memorization problem which was even more difficult
in our case due to the lack of concrete training evaluation methods. We also
did a vast ablation study over many state-of-the-art methods in the
conditional image generation task and developed a stable modified GAN model.

We are still trying to develop methods to overcome the remaining issues. For
example, one idea is to add the prior information of pixel intensity values to
the loss function of the generator like an L2 regularization scheme to
penalize wrong intensity distribution in each batch. In other words, by
minimizing the loss, the model will also minimize the L2 distance between the
generator’s normal output and the generator’s output that has gone through
pixel intensity cut-off. This will be in theory more effective than the hard
cut-off method. We argue that if we only cut off this in a mere post-
processing step, we will lose location-based data in the image for blackening
only the pixel values. Since, despite the incorrectness of pixel values below
7 ADU, they correspond to would-be physical tracks which in principle we need
to have them even with a lower intensity.

Another path worth trying would be to refine the GAN set up to capture the
absolute relation between two layers of the PXD detector like a motion picture
generation using sequence-based GAN models. Furthermore, as we direly need a
reliable training metric, one can create a custom Inception Score (IS) [13],
based on an image classification model trained on simulated events and images
to predict PXD sensor IDs. The original Inception Score is based on the pre-
trained Inception v3 model [29] with which a large number of generated fake
images are classified. Consequently, the probability of the image belonging to
each class is predicted. These predictions are then summarized into the
Inception Score to have a fully automated evaluation metric for the training.
Furthermore, the next step for an ultimate comprehensive validation of
generated hitmaps would be an estimate of a systematic uncertainty on the
tracking efficiency, fake rate, and resolution due using a GAN. If that is
small compared to other systematic uncertainties we can use the model.
Finally, as the PXD background generation is being integrated into the Belle
II software [30], there will be many new experiences that may be very useful
for other experiments as well.

## References

  * [1] T. Abe (Belle II Collaboration), 1011.0352, (2010)
  * [2] K. Akai, K. Furukawa, H. Koiso (SuperKEKB), Nucl. Instrum. Meth. A907, 188 (2018)
  * [3] S. Agostinelli et al. (GEANT4), Nucl.Instrum.Meth. A506, 250 (2003)
  * [4] F. Mueller (DEPFET), JINST 9, C10007 (2014)
  * [5] J. Kemmer, G. Lutz, Nucl. Instrum. Meth. A253, 365 (1987)
  * [6] Goodfellow, Ian J., Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 1406.2661 (2014)
  * [7] de Oliveira, Luke, Michela Paganini, and Benjamin Nachman. Computing and Software for Big Science, 1-24 (2017)
  * [8] Di Sipio, Riccardo, et al. Journal of High Energy Physics 8, 1-17 (2019).
  * [9] Luke Metz, Ben Poole, David Pfau, and Jascha Sohl-Dickstein, International Conference on Learning Representation, pages 1–25 (2017)
  * [10] Jaritz, Maximilian, Raoul De Charette, Emilie Wirbel, Xavier Perrotton, and Fawzi Nashashibi. 2018 International Conference on 3D Vision (3DV), pp. 52-60 (2018)
  * [11] Deng, Jia, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. pp. 248-255 (2009).
  * [12] Brock, Andrew, Jeff Donahue, and Karen Simonyan, 1809.11096 (2018)
  * [13] Salimans, Tim, et al, 1606.03498, (2016)
  * [14] Heusel, Martin, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. arXiv preprint arXiv: 1706.08500 (2017)
  * [15] Zhang, Han, Ian Goodfellow, Dimitris Metaxas, and Augustus Odena. International conference on machine learning, pp. 7354-7363 (2019)
  * [16] Kang, Minguk, and Jaesik Park. 2006.12681 (2020)
  * [17] Vaswani, Ashish, et al. arXiv preprint arXiv:1706.03762 (2017)
  * [18] Zhang, Han, Zizhao Zhang, Augustus Odena, and Honglak Lee. 1910.12027 (2019)
  * [19] Paszke, Adam, et al. arXiv preprint arXiv:1912.01703 (2019)
  * [20] Srebre, Matej, Pascal Schmolz, Hosein Hashemi, Martin Ritter, and Thomas Kuhr. EPJ Web of Conferences 245, p. 02010 (2020)
  * [21] Wu, Yan, Jeff Donahue, David Balduzzi, Karen Simonyan, and Timothy Lillicrap, 1912.00953 (2019)
  * [22] Zhao, Yang, Chunyuan Li, Ping Yu, Jianfeng Gao, and Changyou Chen, 2004.02088 (2020)
  * [23] Miyato, Takeru, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida, 1802.05957 (2018)
  * [24] Saxe, Andrew M., James L. McClelland, and Surya Ganguli. 1312.6120 (2013)
  * [25] Gene Golub and Henk Van der Vorst. Journal of Computational and Applied Mathematics 123, 35–65 (2000)
  * [26] Liu, Kanglin, Wenming Tang, Fei Zhou, and Guoping Qiu. Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 6382-6390 (2019)
  * [27] Mescheder, Lars, Andreas Geiger, and Sebastian Nowozin. International conference on machine learning, pp. 3481-3490 (2018)
  * [28] Agueh, Martial, and Guillaume Carlier. SIAM Journal on Mathematical Analysis 43, no. 2, p904-924 (2011)
  * [29] Szegedy, Christian, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna. Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 2818-2826 (2016)
  * [30] T. Kuhr, C. Pulvermacher, M. Ritter, T. Hauth, N. Braun (Belle-II Framework Software Group), Comput. Softw. Big Sci. 3, 1 (2019)