# Dissipative entanglement generation between two driven qubits in circuit
quantum electrodynamics

Sebastián Luciano Gallardo1, Daniel Domínguez1 and María José Sánchez1,2
1Centro Atómico Bariloche and Instituto Balseiro (Universidad Nacional de
Cuyo), 8400 San Carlos de Bariloche, Río Negro, Argentina.
2Instituto de Nanociencia y Nanotecnología (INN),CONICET-CNEA, Argentina.

###### Abstract

An entangled state generation protocol for a system of two qubits driven with
an ac signal and coupled through a resonator is introduced. We explain the
mechanism of entanglement generation in terms of an interplay between unitary
Landau-Zener-Stückelberg (LZS) transitions induced for appropriate amplitudes
and frequencies of the applied ac signal and dissipative processes dominated
by photon loss. In this way, we found that the steady state of the system can
be tuned to be arbitrarily close to a Bell state, which is independent of the
initial state. Effective two-qubit Hamiltonians that reproduce the resonance
patterns associated with LZS transitions are derived.

††preprint: APS/123-QED

##  I Introduction

The generation and stabilization of entangled states is of fundamental
importance for quantum information applications. In the last two decades,
several proposals have explored strategies based on the use of environmental
noise to obtain and stabilize steady state entanglement Kraus  _et al._
(2008); Verstraete  _et al._ (2009); Tacchino  _et al._ (2018).

Most of these schemes use an external driving field as a tool, with examples
including adiabatic passage protocols Král  _et al._ (2007) \- extensively
employed to generate quantum state transfer Maeda  _et al._ (2006); Zhou  _et
al._ (2017), weak resonant drivings - which enable entanglement stabilization
based on tailoring the relaxation rates in order to generate a highly
entangled steady state Shankar  _et al._ (2013); Leghtas  _et al._ (2013);
Kimchi-Schwartz  _et al._ (2016); Quintana  _et al._ (2013); Campagne-Ibarcq
_et al._ (2018a), or a frequency-modulated signal Li  _et al._ (2020) \- used
to achieve an accelerated formation of dissipative entangled steady states.

These protocols have been tested in several systems such as atomic ensembles
Krauter  _et al._ (2011), trapped ions Barreiro  _et al._ (2011); Lin  _et
al._ (2013); Kienzler  _et al._ (2015), Rydberg atoms Li  _et al._ (2020);
Maeda  _et al._ (2006) and superconducting qubits DiCarlo  _et al._ (2010);
Shankar  _et al._ (2013); Leghtas  _et al._ (2013); Kimchi-Schwartz  _et al._
(2016); Quintana  _et al._ (2013); Campagne-Ibarcq  _et al._ (2018a), to
mention a few.

Recently a mechanism relying on the amplitude-modulation of an ac signal was
proposed to generate steady-state entanglement in a system of two coupled
qubits driven by a large amplitude (non resonant) periodic signal and
interacting with a thermal bath Gramajo  _et al._ (2018, 2021).

Nowadays, circuit quantum electrodynamics (cQED) Blais  _et al._ (2004);
Wallraff  _et al._ (2004); Xiang  _et al._ (2013); Blais  _et al._ (2020,
2021) has been established as one of the leading architectures for studying
quantum computation and quantum simulation, where superconducting qubits are
connected to a transmission line resonator Koch  _et al._ (2007); Gu  _et al._
(2017); Bonifacio  _et al._ (2020). Many important experimental advances have
been achieved in this regard, including the observation of Jaynes-Cummings
ladder Fink  _et al._ (2009) and long-lived qubit-resonator states Paik  _et
al._ (2011), entanglement of distant qubits, realization of one and two qubit
gates and non-demolition readout operations Dewes  _et al._ (2012); DiCarlo
_et al._ (2009, 2010); Stern  _et al._ (2014); Johnson  _et al._ (2010);
Walter  _et al._ (2017); van Loo  _et al._ (2013); Eichler  _et al._ (2012);
Didier  _et al._ (2015); Campagne-Ibarcq  _et al._ (2018b).

In this work we propose a protocol to generate and stabilize maximally
entangled states (in particular, Bell states) in a system of two qubits driven
with an ac signal, which are indirectly coupled via a common resonator.
Although this driving protocol has been implemented in studies of Landau-
Zener-Stückelberg interferometry S.N. Shevchenko (2010); Oliver and Valenzuela
(2009); Ferrón  _et al._ (2012, 2016); Gramajo  _et al._ (2019); Bonifacio
_et al._ (2020) and entanglement generation Gramajo  _et al._ (2018, 2021)
with superconducting qubits, we are not aware of previous proposals employing
ac driven qubits to control entanglement in cQED architectures. Our approach
is rather general and not restricted to the usual weak resonant driving, going
beyond the standard dispersive regime used to couple the resonator for readout
Blais  _et al._ (2020). As it is customary in cQED architectures, we will
assume that the resonator acts as a filter of noise for the qubits Bronn  _et
al._ (2015); Blais  _et al._ (2020), protecting them from spontaneous losses
to the environment. With this in mind, we will model the environment as a
thermal bath coupled to the system mainly through the resonator.

Through an interplay between driving and dissipation, we show that an unique
stationary maximally entangled (Bell) state can be obtained regardless the
initial state of the system, provided the qubits are driven with the
appropriate amplitude and frequency. Moreover, the obtained Bell state is
protected from environmental effects for as long as the driving is applied.

The paper is organized as follows: in Sec.II we do an overview of the system
and present the model Hamiltonian. In Sec.III we solve the unitary driven
dynamics of the system and analize the LZS resonance patterns of the two
relevant transitions involved in the generation of maximally entangled steady
states, once coupling to environment is included in Sec.IV. Additionally two-
qubit Hamiltonians that reproduce the structure of these resonances are also
derived in Sec.III. Conclusions and perspectives are given in Sec.V.

##  II System overview

We study a system composed of two qubits coupled to a bosonic mode within a
resonator, which is itself weakly coupled to a thermal bath with temperature
$T_{\text{b}}$, as is shown schematically in Fig. 1.

Figure 1: Schematic representation of the system under study. The qubit
$i\in\\{1,2\\}$ is periodically driven through $\epsilon_{i}(t)$ and coupled
with strength $g_{i}$ to a mode of frequency $\Omega$ of a resonator. The
resonator is weakly coupled to a thermal bath with strength $\lambda$ via the
bosonic operator $a+a^{\dagger}$.

The Hamiltonian of the qubit $i\in\\{1,2\\}$, including the coupling term to
the resonator is given by

$H_{i}(t)=\frac{\epsilon_{i}(t)}{2}\sigma_{z}^{(i)}+g_{i}(a+a^{\dagger})\sigma_{x}^{(i)},$
(1)

where $\sigma_{j}^{(i)}$ are the Pauli matrices acting on the qubit $i$, and
$a^{\dagger}$ ($a$) is the creation (destruction) operator of the bosonic mode
of the resonator. The qubit $i$’s transition frequency or detuning is
$\epsilon_{i}$ , which can be controlled externally as a function of time. The
coupling of the qubit $i$ to the resonator is of strength $g_{i}$ and we
suppose that the associated operator is transversal to the qubit $i$’s
detuning operator (under this assumption, one can always rotate the qubit
basis such that the coupling to the resonator is through $\sigma_{x}^{(i)}$).
A possible qubit $i$ energy gap induced by a term in Eq.(1) transversal to
$\sigma_{z}^{(i)}$ was neglected under the assumption that it is much smaller
than the corresponding $g_{i}$’s - which is rather justified for several
superconducting qubit systems Oliver  _et al._ (2005). The full Hamiltonian
for the cQED architecture is given by

$\displaystyle H(t)=H_{\text{s}}(t)+H_{\text{b}}+H_{\text{sb}},$ (2)
$\displaystyle H_{\text{s}}(t)=\Omega a^{\dagger}a+\sum_{i=1}^{2}H_{i}(t),$
(3)

where $\Omega$ is the resonator mode frequency. The term $H_{\text{b}}$ in
Eq.(2) represents the bath Hamiltonian, modelled as a continuum of harmonic
oscillators in thermal equilibrium at temperature $T_{\text{b}}$, with ohmic
spectral density $\mathcal{J}(\omega)=\kappa\omega$, where $\kappa$ is a
constant. The term $H_{\text{sb}}$ stands for the interaction between the
system and the thermal bath, which in this work we suppose is through the
operator $(a+a^{\dagger})$ and of strength $\lambda$. The explicit forms of
$H_{\text{b}}$ and $H_{\text{sb}}$ are given in Appendix A.

The driving required for the Bell state generation protocol depends on the
relative sign of the couplings $g_{i}$. We suppose that the couplings are
similar in magnitude but their relative sign could be either equal or
opposite, corresponding to couplings to even or odd modes of the resonator,
respectively. In the following without loss of generality we will consider
couplings with the same sign and the drivings in detuning chosen as

$\epsilon(t)\equiv\epsilon_{1}(t)=\epsilon_{2}(t)=A\cos(\omega t),$ (4)

with $A$ the amplitude and $\omega$ the frequency of the driving. For the case
of opposite coupling signs, the driving should be chosen to be
$\epsilon_{1}(t)=-\epsilon_{2}(t)$. It can be shown that both cases are
related by a local unitary transformation
$H\rightarrow\sigma_{y}^{(2)}H\sigma_{y}^{(2)}$ which keeps the entanglement
generation dynamics invariant.

As we will discuss in detail in Sec.IV, to stimulate the LZS resonances
necessary for this Bell state generation protocol we require that
$0<|\delta_{g}|\ll g_{1}g_{2}/\Omega$, with $\delta_{g}\equiv g_{1}-g_{2}$. As
the relevant involved transitions occur in a timescale $\delta_{g}^{-1}$, the
smaller $\delta_{g}$, the longer it will take to reach the stationary Bell
state.

To solve numerically the dynamics in the purely unitary case (considering only
$H_{\text{s}}(t)$ in Eq.(2)), we diagonalize the evolution operator over a
period of the driving using a 4${}^{\text{th}}$ order Trotter-Suzuki
expansion. In this way we obtain the Floquet states and the associated
quasienergies Shirley (1965). To study the open dynamics we evolved the
system’s density operator using the Floquet-Born-Markov (FBM) master equation
Kohler  _et al._ (1998); Blattmann  _et al._ (2015); Ferrón  _et al._ (2016)
within a moderate Rotating Wave Approximation (RWA), as is detailed in
Appendix B. For the numerical simulations we truncate the Hilbert space to a
finite number of photon levels. We found that retaining the first $5$ photon
levels was sufficient to attain convergence.

##  III Unitary dynamics

In this section we focus on the unitary dynamics described by the Hamiltonian
$H_{\text{s}}(t)$ defined in Eq.(3). It will be useful to define the set of
Bell states of the two qubit system:
$\ket{\Phi_{\pm}}\equiv\frac{1}{\sqrt{2}}(\ket{\uparrow\uparrow}\pm\ket{\downarrow\downarrow})$
and
$\ket{\Psi_{\pm}}\equiv\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}\pm\ket{\downarrow\uparrow})$,
where $\sigma_{z}\ket{\uparrow}=\ket{\uparrow}$, and
$\sigma_{z}\ket{\downarrow}=-\ket{\downarrow}$. The Bell states are maximally
entangled and form a basis for the two qubit Hilbert space.

Figure 2(a) shows the energy spectrum of $H_{\text{s}}$ parametrized as a
function of $\epsilon$, the driving variable. It can be shown that the avoided
crossings at $\epsilon=\pm\Omega$ have energy gaps of first order in $g_{i}$
while for all other integer and half integer values of $\epsilon/\Omega$ the
avoided crossing gaps are of second or higher order in $g_{i}$. Away from all
avoided crossings, the energies and eigenstates of the system satisfy

$\displaystyle H_{\text{s}}(\epsilon)\ket{N\uparrow\uparrow}\ $
$\displaystyle\approx\ (N\Omega+\epsilon)\ket{N\uparrow\uparrow}$ (5)
$\displaystyle H_{\text{s}}(\epsilon)\ket{N\Psi_{\pm}}\ $
$\displaystyle\approx\ \ \ \ \ \ \ \ N\Omega\ket{N\Psi_{\pm}}$ (6)
$\displaystyle H_{\text{s}}(\epsilon)\ket{N\downarrow\downarrow}\ $
$\displaystyle\approx\ (N\Omega-\epsilon)\ket{N\downarrow\downarrow},$ (7)

with $\ket{N}$ the state of the resonator with $N$ photons. In addition, it
can be readily seen from Eq.(3) that for $g_{1}=g_{2}$ ($\delta_{g}=0$), the
singlet states $\ket{N\Psi_{-}}$ are exact eigenstates of
$H_{\text{s}}(\epsilon)$ with energy $N\Omega$. Since they are also
eigenstates of the driving operator
($\propto(\sigma_{z}^{(1)}+\sigma_{z}^{(2)})$), transitions involving the
states $\ket{N\Psi_{-}}$ are forbidden for $\delta_{g}=0$. In fact, defining
$g\equiv(g_{1}+g_{2})/2$ we can write

$\sum_{i}g_{i}\sigma_{x}^{(i)}=g(\sigma_{x}^{(1)}+\sigma_{x}^{(2)})+\frac{\delta_{g}}{2}(\sigma_{x}^{(1)}-\sigma_{x}^{(2)}),$
(8)

and consider the second term in the last equation as a perturbation that
induces transitions between $\ket{N\Psi_{-}}$ and $\ket{(N\pm 1)\Phi_{+}}$,
since we suppose that $\delta_{g}\ll g$. This can be seen from the fact that

$\frac{\delta_{g}}{2}(\sigma_{x}^{(1)}-\sigma_{x}^{(2)})=\delta_{g}\Big{(}\ket{\Psi_{-}}\bra{\Phi_{+}}+\ket{\Phi_{+}}\bra{\Psi_{-}}\Big{)},$
(9)

and that the $(a+a^{\dagger})$ operator either adds or removes a photon from
the resonator state which it acts on.

Figure 2: (a) Energy spectrum of $H_{\text{s}}$ (Eq.(3)) parametrized a
function of $\epsilon/\Omega$, for $g_{1}=0.05\Omega$ and
$g_{2}=0.0485\Omega$. Avoided crossings are located at integer and half-
integer values of $\epsilon/\Omega$. Those at $\epsilon=\pm\Omega$ are of
first order in $g_{i}$ while all others are of second and higher order in
$g_{i}$. Integer values of $\epsilon/\Omega$ are marked with dotted lines. See
text for more details. Energy spectra of the effective Hamiltonians (Eq.(12)
and Eq.(15)) derived to study the transition (b)
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ and the transition (c)
$\ket{0\Psi_{-}}\rightarrow$(other states), respectively. In all cases the
asymptotic eigenstates away from avoided crossings are shown. Figure 3: Upper
panel: Unitary time-averaged transition probabilities calculated numerically
with the Hamiltonian Eq.(3) for the transitions (a)
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ and (c)
$\ket{0\Psi_{-}}\rightarrow$ (other states) , as a function of $A/\omega$ and
$\Omega/\omega$. The patterns show resonances at
$\Omega/\omega=n\in\mathbb{N}$ with width of order $\delta_{g}$ modulated by
the $n^{\text{th}}$ Bessel function $J_{n}(A/\omega)$. Lower panel: Transition
probabilities for (b) $\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ and
(d) $\ket{0\Psi_{-}}\rightarrow$ (other states) , calculated using the
effective two-qubit Hamiltonians Eq.(12) and Eq.(15), respectively. The line
$\Omega/\omega=2$ and the lines for which $A/\omega$ equals the first two
zeros of the $2^{\text{nd}}$ Bessel function are indicated with grey dotted
lines. In all cases the qubit-resonator coupling strengths used are
$g_{1}=0.1\omega$ and $g_{2}=0.097\omega$.

In what follows we consider explicitly transitions involving $\ket{0\Psi_{-}}$
and $\ket{1\Psi_{-}}$ induced by the simultaneous effect of the time dependent
driving $\epsilon(t)$, Eq.(4), and the coupling asymmetry parameter
$\delta_{g}$. These transitions will be relevant for the Bell state generation
mechanism once dissipation is included.

###  III.1 Transitions to $\ket{1\Psi_{-}}$

We begin by studying the transition probability of
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ induced by the periodic
driving. We are interested in this transition because, as we will discuss in
Sec.IV, $\ket{1\Psi_{-}}$ will decay into $\ket{0\Psi_{-}}$ after including
dissipation, which is a maximally entangled state that is stable against
photon loss.

Figure 3(a) shows the time averaged transition probability for
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$, as a function of
$A/\omega$ and $\Omega/\omega$, near a resonance centered around
$\Omega/\omega=2$. Similar resonances of width proportional to $\delta_{g}$
are observed for all integer values of $\Omega/\omega$.

Numerically we find that for these resonances, and for
$\ket{0\uparrow\uparrow}$ as the initial state, the most populated states are
in the subspace $S_{1}$ spanned by
$\\{\ket{0\Phi_{\pm}},\ket{1\Psi_{\pm}}\\}$. Projecting the Hamiltonian of
Eq.(3) into $S_{1}$, the state of the resonator becomes uniquely determined by
the state of the qubits (within $S_{1}$). Thus one can write the terms
involving resonator operators in Eq.(3) in terms of qubit operators:

$\displaystyle a^{\dagger}a|_{S_{1}}=$
$\displaystyle\ket{1\Psi_{+}}\bra{1\Psi_{+}}+\ket{1\Psi_{-}}\bra{1\Psi_{-}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(1-\sigma_{z}^{(1)}\sigma_{z}^{(2)}),$ (10)

and

$\displaystyle\frac{1}{2}(a+a^{\dagger})(\sigma_{x}^{(1)}\pm\sigma_{x}^{(2)})|_{S_{1}}=$
$\displaystyle\ket{0\Phi_{\mp}}\bra{1\Psi_{\pm}}+\ket{1\Psi_{\pm}}\bra{0\Phi_{\mp}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\sigma_{x}^{(1)}\pm\sigma_{x}^{(2)}),$ (11)

where now $\sigma^{(i)}_{j}$ are understood as $4\times 4$ matrices, and
$|_{S_{1}}$ indicates projection into $S_{1}$. Replacing the above expressions
into Eq.(3), one arrives at the effective time dependent Hamiltonian, valid
for studying the transition
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$:

$H_{l}(t)=\sum_{i}\Big{(}\frac{\epsilon(t)}{2}\sigma_{z}^{(i)}+g_{i}\sigma_{x}^{(i)}\Big{)}-\frac{\Omega}{2}\sigma_{z}^{(1)}\sigma_{z}^{(2)}+\frac{\Omega}{2}.$
(12)

Notice that Eq.(12) is the Hamiltonian of two driven qubits coupled
longitudinally studied in Ref.Gramajo  _et al._ , 2017. In the present case
$\Omega$ plays the role of the interaction strength between the qubits and
$2g_{i}$ the role of the intrinsic qubit gaps. Figure 2(b) shows the energy
spectrum of the effective Hamiltonian, Eq.(12), parametrized as a function of
$\epsilon$.

In Fig. 3(b) the time averaged transition probability
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ computed numerically using
Eq.(12) is displayed. The agreement with Fig. 3(a), obtained from the full
Hamiltonian of Eq.(3), is excellent.

For $\delta_{g}^{2}\ll A\omega$ and neglecting the effect of $g$ on the energy
spectrum, the LZS resonance condition for the transition between the two
levels $\ket{0\uparrow\uparrow}$ and $\ket{1\Psi_{-}}$ is $\Omega/\omega=n$
for some integer $n$, where in this case $\Omega$ is playing the same role of
the dc detuning in standard LZS interferometry, as it sets the average energy
difference between the two involved levels. Thus, considering only these two
levels which are separated by a gap of the order of $\delta_{g}$, resonance
patterns centered around $\Omega/\omega=n$ and of width
$\propto\delta_{g}J_{n}(A/\omega)$ (with $J_{n}$ the $n^{\text{th}}$ Bessel
function) are expected Ashhab  _et al._ (2007) and indeed observed in Fig.
3(b). However, the curvature of the resonances is an effect that we found
numerically to be of order $g_{1}g_{2}/\Omega$, as shown in Fig. 4(a). This
effect is not captured under the assumption of neglecting $g$.

Figure 4: Unitary time averaged transition probabilities for the transitions
(a) $\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ and (b)
$\ket{0\Psi_{-}}\rightarrow$(other states) calculated using the effective
Hamiltonians Eq.(12) and Eq.(15) respectively, as a function of $g/\omega$ and
$\Omega/\omega$, for $\delta_{g}=0.003\omega$ and $A=\Omega$. For small $g$,
an approximately quadratic dependence of the resonance displacement as a
function of $g$ is observed, with opposite curvature for both transitions.

###  III.2 Transitions out of $\ket{0\Psi_{-}}$

We now shift our attention to the time averaged probabilities for the
$\ket{0\Psi_{-}}\rightarrow\text{(other states)}$ transition, which is defined
as the sum of all transition probabilities from the initial state
$\ket{0\Psi_{-}}$ to any other states different from it. We are interested in
this transition because it is an indicator of the stability of
$\ket{0\Psi_{-}}$ against unitary transitions induced by the driving that
might take the system out of this state. Figure 3(c) shows numerical results
for the corresponding transition probabilities using the full Hamiltonian
Eq.(3). We have found that the only significantly populated states are in the
subspace $S_{2}$ spanned by
$\\{\ket{0\Psi_{\pm}},\ket{1\Phi_{\pm}},\ket{2\Psi_{+}}\\}$. Thus in the
present case, and unlike the analysis of Sec.III.1, $5$ linearly independent
states are in principle involved in the transition under study.

Numerically we found that the best two-qubit effective model is obtained when
neglecting the state $\ket{2\Psi_{+}}$ in the calculation of the operator
$(a+a^{\dagger})|_{S_{2}}$, and the state $\ket{0\Psi_{+}}$ in the calculation
of the operator $a^{\dagger}a|_{S_{2}}$. Doing this approximation, one obtains

$\displaystyle a^{\dagger}a|_{S_{2}}\approx$
$\displaystyle\ket{1\Phi_{+}}\bra{1\Phi_{+}}+\ket{1\Phi_{-}}\bra{1\Phi_{-}}+2\ket{2\Psi_{+}}\bra{2\Psi_{+}}$
$\displaystyle=$ $\displaystyle
1+\frac{1}{2}(\sigma_{x}^{(1)}\sigma_{x}^{(2)}+\sigma_{y}^{(1)}\sigma_{y}^{(2)}),$
(13)

$\displaystyle\frac{1}{2}(a+a^{\dagger})(\sigma_{x}^{(1)}\pm\sigma_{x}^{(2)})|_{S_{2}}\approx$
$\displaystyle\ket{0\Psi_{\pm}}\bra{1\Phi_{\mp}}+\ket{1\Phi_{\mp}}\bra{0\Psi_{\pm}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\sigma_{x}^{(1)}\pm\sigma_{x}^{(2)}),$ (14)

and we arrive at the effective Hamiltonian, valid for studying the transition
$\ket{0\Psi_{-}}\rightarrow$(other states):

$\displaystyle H_{\text{tr}}(t)=$
$\displaystyle\sum_{i}\Big{(}\frac{\epsilon(t)}{2}\sigma_{z}^{(i)}+g_{i}\sigma_{x}^{(i)}\Big{)}+$
$\displaystyle\frac{\Omega}{2}(\sigma_{x}^{(1)}\sigma_{x}^{(2)}+\sigma_{y}^{(1)}\sigma_{y}^{(2)})+\Omega,$
(15)

which is the Hamiltonian of two transversally coupled and symmetrically driven
qubits, with $\Omega$ playing the role of the interaction strength and
$2g_{i}$ the role of the intrinsic qubit gaps Gramajo  _et al._ (2017). Figure
2(c) shows its energy spectrum parametrized as a function of $\epsilon$.
Numerical results for the transition probability
$\ket{0\Psi_{-}}\rightarrow\text{(other states)}$ computed from Eq.(15) are
shown in Fig.3(d). Except for the additional flatter resonances around integer
$\Omega/\omega$, which correspond to resonances to states outside $S_{2}$, the
agreement with Fig.3(c), obtained using the full Hamiltonian Eq.(3) is
remarkable.

The LZS resonance condition for both transitions,
$\ket{0\Psi_{-}}\rightarrow\ket{1\uparrow\uparrow}$ and
$\ket{0\Psi_{-}}\rightarrow\ket{1\downarrow\downarrow}$ respectively, is again
$\Omega/\omega=n$ for some integer $n$. Resonances around these values of
$\Omega/\omega$ of width $\propto\delta_{g}J_{n}(A/\omega)$ are expected and
observed in Fig. 3(d), in analogy to the results of Sec.III.1. The effect of
$g$ on the curvature of the resonances out of $\ket{0\Psi_{-}}$ is again of
order $g_{1}g_{2}/\Omega$ but of opposite sign to that of the transitions to
$\ket{1\Psi_{-}}$, as can be seen in Fig. 4(b). This implies that there are
values of the driving parameters, $A$ and $\omega$, for which the transitions
to $\ket{1\Psi_{-}}$ are stimulated but those out of $\ket{0\Psi_{-}}$ are
not. This is the key point that will be made use of to generate Bell states
once dissipation is included, as is explained below.

##  IV Dissipation induced Bell state generation

In the previous section we concluded that starting from the initial state
$\ket{0\uparrow\uparrow}$, there are regions in the $A/\omega-\Omega/\omega$
plane where a unitary resonance to $\ket{1\Psi_{-}}$ is stimulated but no
unitary resonance involving $\ket{0\Psi_{-}}$ is so. When the driving
amplitude A and frequency $\omega$ are chosen as to select one of these
points, the process

$\ket{0\uparrow\uparrow}\stackrel{{\scriptstyle
H_{\text{s}}(t)}}{{\rightarrow}}\ket{1\Psi_{-}}\stackrel{{\scriptstyle\text{PL}}}{{\rightarrow}}\ket{0\Psi_{-}}$
(16)

can occur once dissipation is included, where the first transition is unitary
and induced by the driven Hamiltonian $H_{s}(t)$, and the second one is the
loss of a photon to the environment. The same process can also take place when
starting from $\ket{0\downarrow\downarrow}$ or $\ket{0\Psi_{+}}$, since they
show similar unitary resonances for
$\ket{0\downarrow\downarrow}\rightarrow\ket{1\Psi_{-}}$ and
$\ket{0\Psi_{+}}\rightarrow\ket{1\Psi_{-}}$.

Since the open system dynamics is linear due to the assumed weak coupling to
the environment, one can understand the evolution of the system’s density
matrix as the independent evolution of its ensemble members which, for
$T_{\text{b}}\ll\Omega$, will eventually reach one of the states with zero
photons. If the reached states are different from $\ket{0\Psi_{-}}$, they will
go through the process defined in Eq.(16) ending up at least partially in
$\ket{0\Psi_{-}}$. As all transitions involving $\ket{0\Psi_{-}}$ are out of
resonance, gradually the population of the system’s density matrix accumulates
in this state.

Figure 5: Time averaged (a) population of $\ket{0\Psi_{-}}$ and (b)
concurrence of the steady state of the driven dissipative system of two qubits
coupled to a resonator, obtained after solving numerically the FBM master
equation associated to the Hamiltonian defined in Eq.(2). Asymmetrical
resonance patterns are observed at integer values of $\Omega/\omega$. The
coupling strengths are the same as in Fig. 3, with bath parameters
$\kappa\lambda^{2}=0.0001$ and $T_{\text{b}}=0.001\omega$. The line
$\Omega/\omega=2$ and the lines for which $A/\omega$ equals the first two
zeros of the $2^{\text{nd}}$ Bessel function are indicated with grey dotted
lines.

Figure 5(a) shows numerical results for the time averaged population of
$\ket{0\Psi_{-}}$ in the stationary state as a function of $A/\omega$ and
$\Omega/\omega$ obtained after solving numerically the FBM master equation for
the system’s density matrix $\rho$ (see Appendix B for details). The
stationary state is found to be unique and $T$-periodic in time. The behavior
of the stationary population closely follows what is predicted by the previous
argument: population maxima very close to $1$ are observed for points in the
parameter space along the unitary resonance patterns of
$\ket{0\uparrow\uparrow}\rightarrow\ket{1\Psi_{-}}$ and population minima,
very close to $0$, are obtained for points along the unitary
$\ket{0\Psi_{-}}\rightarrow\text{(other states)}$ resonances.

To quantify the degree of entanglement of the qubits, we use the concurrence
as a measure Wootters (1998),

$\mathcal{C}[\rho_{\text{q}}]=\max\\{0,r_{3}-r_{2}-r_{1}-r_{0}\\},$ (17)

defined in terms of the qubits’ density matrix
$\rho_{\text{q}}=\text{Tr}_{\text{r}}(\rho)$, where the trace operation is
over the states of the resonator and $r_{i}$ are the real-valued eigenvalues
of

$\sqrt{\sqrt{\rho_{\text{q}}}\sigma_{y}^{(1)}\sigma_{y}^{(2)}\rho_{\text{q}}^{*}\sigma_{y}^{(1)}\sigma_{y}^{(2)}\sqrt{\rho_{\text{q}}}},$
(18)

sorted in ascending order, where $\rho_{\text{q}}^{*}$ is the complex
conjugate of $\rho_{\text{q}}$ and the conjugation must be done in a separable
basis. The concurrence takes a value of $0$ for a separable state, a value of
$1$ for a maximally entangled state, and values in between for partially
entangled states.

Figure 5(b) shows the time averaged concurrence of the stationary state. It is
observed that the maxima of the concurrence that are close to $1$ are achieved
only when the state $\ket{0\Psi_{-}}$ is populated, indicating that this is
the only maximally entangled state that is generated. It is also noteworthy
that for driving parameters lying outside the mentioned resonances, which
constitute the vast majority of the points in the $A/\omega-\Omega/\omega$
plane (including the case of no driving at all $A=0$), the stationary state is
either separable or almost separable.

As we have already mentioned, to attain $\ket{0\Psi_{-}}\bra{0\Psi_{-}}$ as
the stationary state of the system, a necessary condition is to find points in
the plane $A/\omega-\Omega/\omega$ where the resonance conditions for the
transitions to $\ket{1\Psi_{-}}$ and out of $\ket{0\Psi_{-}}$ do not overlap.
This can only be fulfilled if the maximum resonance deviation from the
condition $\Omega/\omega=n$ (of order $g_{1}g_{2}/\Omega$) is much greater
than the resonance width (of order $\delta_{g}$), and this the source of the
requirement $\delta_{g}\ll g_{1}g_{2}/\Omega$ . The optimal entanglement
generation (maximal area and intensity of concurrence patterns) is achieved
for amplitudes $A=\Omega$ and frequencies $\omega$ such that $\Omega$ is below
an integer multiple of $\omega$ by a frequency of the order of
$g_{1}g_{2}/\Omega$, i.e. $A=\Omega=n\omega-\mathcal{O}(g_{1}g_{2}/\Omega)$.

Finally, Fig.6 shows the temporal dynamics of relevant populations of the
system’s density matrix at a point of high entanglement generation
$A=\Omega=1.983\omega$, for the system starting in the state
$\rho_{0}=\ket{0\uparrow\uparrow}\bra{0\uparrow\uparrow}$. It is seen that
even though the dynamics of the populations is complicated, a clear resonance
to the state $\ket{1\Psi_{-}}$ is stimulated in a timescale of the order of
$\delta_{g}^{-1}$, which gradually decays into $\ket{0\Psi_{-}}$ via photon
loss to the environment. Population accumulates in this state, and given
enough time the system ends up essentially at
$\rho_{\infty}=\ket{0\Psi_{-}}\bra{0\Psi_{-}}$.

We thus conclude that by applying a driving with appropriate amplitude and
frequency it is possible to populate the maximally entangled state
$\ket{0\Psi_{-}}$ independently of the initial state of the system.

Figure 6: Temporal dependence of the system’s density matrix populations,
starting from the initial state
$\ket{0\uparrow\uparrow}\bra{0\uparrow\uparrow}$, for up to (a) $1000$ (linear
scale) and (b) $10^{6}$ (log scale) driving periods, dissipation included. The
driving parameters are $A=\Omega=1.983\omega$, and the rest of the parameters
are the same as in Fig.5. The populations of $\ket{0\uparrow\uparrow}$,
$\ket{0\downarrow\downarrow}$, $\ket{0\Psi_{-}}$ and $\ket{1\Psi_{-}}$ are
shown. It is seen that for short times a resonance between
$\ket{0\uparrow\uparrow}$, $\ket{0\downarrow\downarrow}$ and $\ket{1\Psi_{-}}$
is stimulated. For later times and via photon loss population gradually
accumulates in $\ket{0\Psi_{-}}$.

##  V Conclusions

In this work we have presented an entanglement generation protocol for a
system of two qubits coupled through a resonator. A maximally entangled steady
state is achieved when a symmetric ac driving is applied over both qubits with
appropriate amplitude and frequency.

The areas of entanglement generation in the $A/\omega$-$\Omega/\omega$ plane
are associated to resonant unitary transitions into $\ket{1\Psi_{-}}$ and out
of $\ket{0\Psi_{-}}$. If the driving is such that transitions to
$\ket{1\Psi_{-}}$ are in resonance, but transitions out of $\ket{0\Psi_{-}}$
are not, once dissipation is included the system will accumulate population in
$\ket{0\Psi_{-}}$ via photon loss from $\ket{1\Psi_{-}}$. All the relevant
features of the unitary resonance patterns were described in terms of
effective Hamiltonians for two driven qubits.

The optimal steady state entanglement generation (in terms of maximal area and
intensity) is attained for driving amplitudes $A=\Omega$ and frequencies such
that $\Omega/\omega$ is slightly below an integer number.

The proposed protocol shows advantages with respect to other methods for
entanglement generation. Due to dissipation, the steady state is independent
of the initial state of the system Kraus  _et al._ (2008); Verstraete  _et
al._ (2009) and, once the entangled state is obtained, it is protected from
environmental effects for as long as the driving is applied. Moreover, the
present proposal allows entangling distant and strongly driven qubits which
are, for example, a microwave waveguide apart. Therefore, it is expected that
our scheme could add a robust means to realize entanglement protocols in
setups extensively used nowadays in cQED Blais  _et al._ (2020, 2021).

## Acknowlegments

We acknowledge support from CNEA, CONICET (PIP11220150100756), ANPCyT
(PICT2016-0791 and PICT2019-0654) and UNCuyo (06/C591).

##  Appendix A Model of the environment

We model the bath and its interaction with the system using the Caldeira-
Leggett model Breuer and Petruccione (2006):

$\displaystyle H_{\text{b}}=$ $\displaystyle\int_{0}^{\infty}d\omega\omega
b_{\omega}^{\dagger}b_{\omega}$ (19) $\displaystyle H_{\text{sb}}=$
$\displaystyle\lambda
A\int_{0}^{\infty}d\omega\sqrt{\mathcal{J}(\omega)}(b_{\omega}+b_{\omega}^{\dagger})+H_{\text{rn}}$
(20) $\displaystyle H_{\text{rn}}=$
$\displaystyle\lambda^{2}A^{2}\int_{0}^{\infty}d\omega\frac{\mathcal{J}(\omega)}{\omega},$
(21)

where $b_{\omega}$ and $b_{\omega}^{\dagger}$ are the creation and destruction
operators of the harmonic oscillator continuum, $A$ is a unitless system
operator, $\lambda$ is a coupling strength, and $H_{\text{rn}}$ is a
renormalization term to cancel all Lamb shifts induced on the system by the
thermal bath.

For the system under study we chose $A=a+a^{\dagger}$ and
$\mathcal{J}(\omega)=\kappa\omega$, with $\kappa$ a constant with units of
energy-2. Usually, it is necessary to put a cutoff frequency in the bath
spectral density, but for our purposes, and since we have already cancelled
out the Lamb shifts, it is permissible to take this cutoff frequency as
infinity (greater than all energy scales of the problem) as we have implicitly
done by the choice of $\mathcal{J}(\omega)$.

##  Appendix B Floquet-Born-Markov master equation

Floquet theory is widely used to study time periodic unitary quantum systems
Shirley (1965). It shows that for a quantum system with a $T$-periodic
Hamiltonian $H(t)$, all solutions are a linear combination of a single basis
of states of the form $e^{-i\epsilon_{\alpha}t}\ket{u_{\alpha}(t)}$. Here,
$\ket{u_{\alpha}(t)}$ is $T$-periodic and is called a Floquet state, and
$\epsilon_{\alpha}$ is called its corresponding quasienergy. To find the
Floquet states and their quasienergies it is customary to diagonalize the
evolution operator $U(t,t_{0})$ over a period of the driving, since the
Floquet states satisfy the eigenvalue equation:

$U(t+T,t)\ket{u_{\alpha}(t)}=e^{-i\epsilon_{\alpha}T}\ket{u_{\alpha}(t)}.$
(22)

The Floquet-Born-Markov master equation Kohler  _et al._ (1998); Ferrón  _et
al._ (2012); Blattmann  _et al._ (2015) allows modelling dissipative processes
in periodically driven systems for sufficiently weak coupling to the
environment. It is a linear Markovian differential equation for the time
evolution of the system’s density matrix:

$\partial_{t}\rho_{\alpha\beta}(t)=-i(\epsilon_{\alpha}-\epsilon_{\beta})\rho_{\alpha\beta}(t)+\sum_{\alpha^{\prime}\beta^{\prime}}L_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(t)\rho_{\alpha^{\prime}\beta^{\prime}}(t),$
(23)

where $\rho_{\alpha\beta}(t)=\bra{u_{\alpha}(t)}\rho(t)\ket{u_{\beta}(t)}$ are
the components of the system’s density matrix in a Floquet basis. The first
term in (23) corresponds to the unitary evolution of the system, while the
second one takes into account dissipative effects. The transition rates
$L_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(t)$ are $T$-periodic and can be
Fourier expanded as

$L_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(t)=\sum_{q}L^{q}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}e^{-iq\omega
t},$ (24)

with $q\in\mathbb{Z}$, and where the coefficients
$L^{q}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}$ are given by

$\displaystyle
L^{q}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}=\lambda^{2}\sum_{k}\Big{(}g_{\alpha\alpha^{\prime}}^{k}A_{\alpha\alpha^{\prime}}^{k}A_{\beta^{\prime}\beta}^{-k-q}+g_{\beta\beta^{\prime}}^{k}A_{\alpha\alpha^{\prime}}^{k-q}A_{\beta^{\prime}\beta}^{-k}$
$\displaystyle-\sum_{\eta}(\delta_{\beta\beta^{\prime}}g_{\eta\alpha^{\prime}}^{k}A_{\alpha\eta}^{-k-q}A_{\eta\alpha^{\prime}}^{k}+\delta_{\alpha\alpha^{\prime}}g_{\eta\beta^{\prime}}^{k}A_{\beta^{\prime}\eta}^{-k}A_{\eta\beta}^{k-q})\Big{)}.$
(25)

In the last expression, the index $\eta$ runs over the indices of the Floquet
basis, $\delta_{\alpha\beta}$ is the Kronecker delta and we defined

$A_{\alpha\beta}^{q}=\sum_{k}\bra{u_{\alpha}^{k}}A\ket{u_{\beta}^{q+k}}$ (26)

and

$g_{\alpha\beta}^{k}=g(\epsilon_{\alpha}-\epsilon_{\beta}+k\omega).$ (27)

In Eq.(26), $\ket{u_{\alpha}^{k}}$ is the $k$th Fourier component of the
Floquet state $\ket{u_{\alpha}(t)}$. In Eq. (27), $g(\omega^{\prime})$ is the
Fourier transformed correlation function of the thermal bath, which can be
expressed in terms of its spectral density and the Bose occupation number
$n_{T_{\text{b}}}(\omega^{\prime})=1/(e^{\omega^{\prime}/T_{\text{b}}}-1)$ as

$g(\omega^{\prime})=\begin{cases}\mathcal{J}(\omega^{\prime})n_{T_{\text{b}}}(\omega^{\prime})\
\ \ &\Leftarrow\omega^{\prime}>0\\\
-\mathcal{J}(-\omega^{\prime})n_{T_{\text{b}}}(\omega^{\prime})&\Leftarrow\omega^{\prime}<0,\end{cases}$
(28)

with the value of $g$ at $\omega^{\prime}=0$ obtained by taking the
appropriate limit.

For sufficiently weak coupling to the environment, such that the maximum rate
of relaxation or decoherence is much smaller than the driving frequency, a
(moderate) RWA is justified in the transition rates, Eq.(24). This sets the
terms with $q\neq 0$ effectively to zero, yielding the simplified expression:

$L_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\approx
R_{\alpha\beta\alpha^{\prime}\beta\beta^{\prime}}+R_{\beta\alpha\beta^{\prime}\alpha^{\prime}}^{*}-\sum_{\eta}(\delta_{\beta\beta^{\prime}}R_{\eta\eta\alpha^{\prime}\alpha}+\delta_{\alpha\alpha^{\prime}}R_{\eta\eta\beta^{\prime}\beta}^{*}),$
(29)

in terms of the quantities

$R_{\alpha\beta\alpha^{\prime}\beta^{\prime}}=\sum_{k}g_{\alpha\alpha^{\prime}}^{k}A_{\alpha\alpha^{\prime}}^{k}A_{\beta\beta^{\prime}}^{-k}.$
(30)

With this approximation, Eq.(23) no longer depends explicitly on time in the
Floquet basis.

For the cases studied in this work, it was found that the operator
$\Lambda_{\alpha\beta\alpha^{\prime}\beta^{\prime}}=-i(\epsilon_{\alpha}-\epsilon_{\beta})\delta_{\alpha\alpha^{\prime}}\delta_{\beta\beta^{\prime}}+L_{\alpha\beta\alpha^{\prime}\beta^{\prime}}$
can be numerically diagonalized in terms of left and right eigenvectors which
are density matrices. That is,

$\Lambda\rho_{\mu}^{\text{R}}=\rho_{\mu}^{\text{L}}\Lambda=\zeta_{\mu}\rho_{\mu},$
(31)

with $\zeta_{\mu}\in\mathbb{C}$ and
$\frac{1}{N}\text{Tr}(\rho_{\mu}^{\text{L}}\rho_{\nu}^{\text{R}})=\delta_{\mu\nu}$,
where $N$ is the dimension of the system’s Hilbert space. Once these
eigenvectors and eigenvalues are obtained, density matrices can be evolved
readily by projecting on this eigensystem,

$\rho(t)=\sum_{\mu}c_{\mu}e^{\zeta_{\mu}(t-t_{0})}\rho_{\mu}^{\text{R}},\ \
c_{\mu}=\frac{1}{N}\text{Tr}(\rho_{\mu}^{\text{L}}\rho(t_{0})).$ (32)

The real parts of $\zeta_{\mu}$ (which are always negative) are the
decoherence and relaxation rates. The stationary state $\rho_{\infty}$, which
for all cases studied in this work can be found and is unique, is defined (in
the Floquet basis) as the state $\rho^{\text{R}}_{\mu}$ with $\zeta_{\mu}=0$.
It is constant in the Floquet basis and therefore $T$-periodic in the original
system basis.

To calculate time averaged functions $f(\rho)$ of the system’s density matrix
in the stationary state, such as populations or concurrence, we make use of
the periodicity of $\rho_{\infty}$ and numerically integrate

$\bar{f}=\frac{1}{T}\int_{0}^{T}dtf(\rho_{\infty}(t)).$ (33)

## References

  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. Ignacio Cirac,  5, 633 EP (2009).
  * Tacchino  _et al._ (2018) F. Tacchino, A. Auffèves, M. F. Santos, and D. Gerace, Phys. Rev. Lett. 120, 063604 (2018).
  * Král  _et al._ (2007) P. Král, I. Thanopulos, and M. Shapiro, Rev. Mod. Phys. 79, 53 (2007).
  * Maeda  _et al._ (2006) H. Maeda, J. H. Gurian, D. V. L. Norum, and T. F. Gallagher, Phys. Rev. Lett. 96, 073002 (2006).
  * Zhou  _et al._ (2017) B. B. Zhou, A. Baksic, H. Ribeiro, C. Yale, F. Heremans, P. Jerger, A. Auer, G. Burkard, A. Clerk, and D. Awschalom, Nature Physics 13, 330 (2017).
  * Shankar  _et al._ (2013) S. Shankar, M. Hatridge, Z. Leghtas, K. Sliwa, A. Narla, U. Vool, S. M. Girvin, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Nature 504, 419 (2013).
  * Leghtas  _et al._ (2013) Z. Leghtas, U. Vool, S. Shankar, M. Hatridge, S. M. Girvin, M. H. Devoret, and M. Mirrahimi, Physical Review A 88, 023849 (2013).
  * Kimchi-Schwartz  _et al._ (2016) M. E. Kimchi-Schwartz, L. Martin, E. Flurin, C. Aron, M. Kulkarni, H. E. Tureci, and I. Siddiqi, Phys. Rev. Lett. 116, 240503 (2016).
  * Quintana  _et al._ (2013) C. M. Quintana, K. D. Petersson, L. W. McFaul, S. J. Srinivasan, A. A. Houck, and J. R. Petta, Phys. Rev. Lett. 110, 173603 (2013).
  * Campagne-Ibarcq  _et al._ (2018a) P. Campagne-Ibarcq, E. Zalys-Geller, A. Narla, S. Shankar, P. Reinhold, L. Burkhart, C. Axline, W. Pfaff, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Phys. Rev. Lett. 120, 200501 (2018a).
  * Li  _et al._ (2020) R. Li, D. Yu, S.-L. Su, and J. Qian, Phys. Rev. A 101, 042328 (2020).
  * Krauter  _et al._ (2011) H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, Phys. Rev. Lett. 107, 080503 (2011).
  * Barreiro  _et al._ (2011) J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 EP (2011).
  * Lin  _et al._ (2013) Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler, A. S. Sørensen, D. Leibfried, and D. J. Wineland,  504, 415 EP (2013).
  * Kienzler  _et al._ (2015) D. Kienzler, H.-Y. Lo, B. Keitch, L. de Clercq, F. Leupold, F. Lindenfelser, M. Marinelli, V. Negnevitsky, and J. P. Home, Science 347, 53 (2015).
  * DiCarlo  _et al._ (2010) L. DiCarlo, M. D. Reed, L. Sun, B. R. Johnson, J. M. Chow, J. M. Gambetta, L. Frunzio, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Nature 467, 574 (2010).
  * Gramajo  _et al._ (2018) A. L. Gramajo, D. Domínguez, and M. J. Sánchez, Phys. Rev. A 98, 042337 (2018).
  * Gramajo  _et al._ (2021) A. L. Gramajo, D. Domínguez, and M. J. Sánchez, Phys. Rev. A 104, 032410 (2021).
  * Blais  _et al._ (2004) A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 69, 062320 (2004).
  * Wallraff  _et al._ (2004) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nature 431, 162 (2004).
  * Xiang  _et al._ (2013) Z.-L. Xiang, S. Ashhab, J. You, and F. Nori, Reviews of Modern Physics 85, 623 (2013).
  * Blais  _et al._ (2020) A. Blais, S. M. Girvin, and W. D. Oliver, Nature Physics 16, 247 (2020).
  * Blais  _et al._ (2021) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Koch  _et al._ (2007) J. Koch, T. Yu, M. Terri, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. Devoret, S. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * Gu  _et al._ (2017) X. Gu, A. F. Kockum, A. Miranowicz, Y. xi Liu, and F. Nori, Physics Reports 718-719, 1 (2017), microwave photonics with superconducting quantum circuits.
  * Bonifacio  _et al._ (2020) M. Bonifacio, D. Domínguez, and M. J. Sánchez, Physical Review B 101, 245415 (2020).
  * Fink  _et al._ (2009) J. M. Fink, R. Bianchetti, M. Baur, M. Göppl, L. Steffen, S. Filipp, P. J. Leek, A. Blais, and A. Wallraff, Phys. Rev. Lett. 103, 083601 (2009).
  * Paik  _et al._ (2011) H. Paik, D. I. Schuster, L. S. Bishop, G. Kirchmair, G. Catelani, A. P. Sears, B. R. Johnson, M. J. Reagor, L. Frunzio, L. I. Glazman, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Phys. Rev. Lett. 107, 240501 (2011).
  * Dewes  _et al._ (2012) A. Dewes, R. Lauro, F. R. Ong, V. Schmitt, P. Milman, P. Bertet, D. Vion, and D. Esteve, Phys. Rev. B 85, 140503 (2012).
  * DiCarlo  _et al._ (2009) L. DiCarlo, J. M. Chow, J. M. Gambetta, L. S. Bishop, B. R. Johnson, D. Schuster, J. Majer, A. Blais, L. Frunzio, S. Girvin, _et al._ , Nature 460, 240 (2009).
  * Stern  _et al._ (2014) M. Stern, G. Catelani, Y. Kubo, C. Grezes, A. Bienfait, D. Vion, D. Esteve, and P. Bertet, Phys. Rev. Lett. 113, 123601 (2014).
  * Johnson  _et al._ (2010) B. R. Johnson, M. D. Reed, A. A. Houck, D. I. Schuster, L. S. Bishop, E. Ginossar, J. M. Gambetta, L. DiCarlo, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, Nature Physics 6, 663 (2010).
  * Walter  _et al._ (2017) T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, Phys. Rev. Applied 7, 054020 (2017).
  * van Loo  _et al._ (2013) A. F. van Loo, A. Fedorov, K. Lalumiere, B. C. Sanders, A. Blais, and A. Wallraff, Science 342, 1494–1496 (2013).
  * Eichler  _et al._ (2012) C. Eichler, C. Lang, J. M. Fink, J. Govenius, S. Filipp, and A. Wallraff, Physical Review Letters 109 (2012), 10.1103/physrevlett.109.240501.
  * Didier  _et al._ (2015) N. Didier, J. Bourassa, and A. Blais, Physical Review Letters 115 (2015), 10.1103/physrevlett.115.203601.
  * Campagne-Ibarcq  _et al._ (2018b) P. Campagne-Ibarcq, E. Zalys-Geller, A. Narla, S. Shankar, P. Reinhold, L. Burkhart, C. Axline, W. Pfaff, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Phys. Rev. Lett. 120, 200501 (2018b).
  * S.N. Shevchenko (2010) F. N. S.N. Shevchenko, S. Ashhab, Physics Reports 492, 1 (2010).
  * Oliver and Valenzuela (2009) W. D. Oliver and S. O. Valenzuela, Quantum Information Processing 8, 261 (2009).
  * Ferrón  _et al._ (2012) A. Ferrón, D. Domínguez, and M. J. Sánchez, Phys. Rev. Lett. 109, 237005 (2012).
  * Ferrón  _et al._ (2016) A. Ferrón, D. Domínguez, and M. J. Sánchez, Phys. Rev. B 93, 064521 (2016).
  * Gramajo  _et al._ (2019) A. L. Gramajo, D. Domínguez, and M. J. Sánchez, Phys. Rev. B 100, 075410 (2019).
  * Bronn  _et al._ (2015) N. T. Bronn, E. Magesan, N. A. Masluk, J. M. Chow, J. M. Gambetta, and M. Steffen, IEEE Transactions on Applied Superconductivity 25, 1 (2015).
  * Oliver  _et al._ (2005) W. D. Oliver, Y. Yu, J. C. Lee, K. K. Berggren, L. S. Levitov, and T. P. Orlando, Science 310, 1653 (2005).
  * Shirley (1965) J. H. Shirley, Physical Review 138, B979 (1965).
  * Kohler  _et al._ (1998) S. Kohler, R. Utermann, P. Hänggi, and T. Dittrich, Phys. Rev. E 58, 7219 (1998).
  * Blattmann  _et al._ (2015) R. Blattmann, P. Hänggi, and S. Kohler, Phys. Rev. A 91, 042109 (2015).
  * Gramajo  _et al._ (2017) A. L. Gramajo, D. Domínguez, and M. J. Sánchez, Eur. Phys. J. B 90, 255 (2017).
  * Ashhab  _et al._ (2007) S. Ashhab, J. R. Johansson, A. M. Zagoskin, and F. Nori, Phys. Rev. A 75, 063414 (2007).
  * Wootters (1998) W. K. Wootters, Physical Review Letters 80, 2245 (1998).
  * Breuer and Petruccione (2006) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2006).


# Connecting anomaly and tunneling methods for Hawking effect through
chirality

Rabin Banerjee, Bibhas Ranjan Majhi
S. N. Bose National Centre for Basic Sciences,
JD Block, Sector III, Salt Lake, Kolkata-700098, India
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The role of chirality is discussed in unifying the anomaly and the tunneling
formalisms for deriving the Hawking effect. Using the chirality condition and
starting from the familiar form of the trace anomaly, the chiral
(gravitational) anomaly, manifested as a nonconservation of the stress tensor,
near the horizon of a black hole, is derived. Solution of this equation yields
the stress tensor whose asymptotic infinity limit gives the Hawking flux.
Finally, use of the same chirality condition in the tunneling formalism gives
the Hawking temperature that is compatible with the flux obtained by anomaly
method.

Introduction: Ever since Hawking’s original observation [1] that black holes
radiate, there have been several derivations [2, 3, 4, 5, 6, 7, 8] of this
effect. A common feature in these derivations is the universality of the
phenomenon; the Hawking radiation is determined universally by the horizon
properties of the black hole leading to the same answer. This, in the absence
of direct experimental evidence, definitely reinforces Hawking’s original
conclusion. Moreover, it strongly suggests that there is some fundamental
mechanism which could, in some sense, unify the various approaches.

In this paper we show that chirality is the common property which connects the
tunneling formalism [4, 5, 9, 10] and the anomaly method [3, 6, 7, 8, 11, 12,
13] in studying Hawking effect. Apart from being among the most widely used
approaches, interest in both the anomaly and tunneling methods has been
revived recently leading to different variations and refinements in them [7,
8, 9, 10, 11, 12, 13, 14]. The calculation will be performed using a family of
metrics that includes a subset of the stationary, spherically symmetric space-
times which are asymptotically flat. Also, the results are derived using
mostly physical reasoning and do not require any specific technical skill.

Before commencing on our analysis we briefly recapitulate the basic tenets of
the tunneling and anomaly methods. The idea of a tunneling description, quite
akin to what we know in usual quantum mechanics where classically forbidden
processes might be allowed through quantum tunneling, dates back to 1976 [15].
Present day computations generally follow either the null geodesic method [5,
9] or the Hamilton-Jacobi method [4, 9, 10], both of which rely on the
semiclassical WKB approximation yielding equivalent results. The essential
idea is that a particle-antiparticle pair forms close to the event horizon.
The ingoing negative energy mode is trapped inside the horizon while the
outgoing positive energy mode is observed at infinity as the Hawking flux.

Although the notion of an anomaly, which represents the breakdown of some
classical symmetry upon quantisation, is quite old, its implications for
Hawking effect were first studied in [3]. It was based on the conformal
(trace) anomaly but the findings were confined only to two dimensions. However
it is possible to apply this method to general dimensions. Recently a new
method was put forward in [6] where a general (any dimensions) derivation was
given. It was based on the well known fact that the effective theory near the
event horizon is a two dimensional conformal theory. The ingoing modes are
trapped within the horizon and cannot contribute to the effective theory near
the horizon. Thus the near horizon theory becomes a two dimensional chiral
theory. Such a chiral theory suffers from a general coordinate
(diffeomorphism) anomaly manifested by a nonconservation of the stress tensor.
Using this gravitational anomaly and a suitable boundary condition the Hawking
flux was obtained. A simplified version of this method was given in [7]. This
was followed by another, new, anomaly based approach in [8, 11].

The first step in our procedure is to derive the gravitational anomaly using
the notion of chirality. This is a new method of obtaining the gravitational
anomaly. Once this anomaly is obtained, the flux is easily deduced. Exploiting
the same notion of chirality the probability of the outgoing mode in the
tunneling approach will be computed. The Hawking temperature then follows from
this probability. At an intermediate stage of this computation we further show
that the chiral modes obtained in the tunneling formalism reproduce the
gravitational anomaly thereby completing the circle of arguments regarding the
connection of the two approaches.

Metric and null coordinates: Consider a black hole characterised by a
spherically symmetric, static space-time and asymptotically flat metric of the
form,

$\displaystyle ds^{2}=F(r)dt^{2}-\frac{dr^{2}}{F(r)}-r^{2}d\Omega^{2}$ (1)

whose event horizon $r=r_{H}$ is defined by $F(r_{H})=0$. Now it is well known
[16, 6] that near the event horizon the effective theory reduces to a two
dimensional conformal theory whose metric is given by the ($r-t$) sector of
the original metric (1).

It is convenient to express (1) in the null tortoise coordinates which are
defined as,

$\displaystyle u=t-r_{*},\,\,\,v=t+r_{*};\,\,\,\,dr_{*}=\frac{dr}{F(r)}.$ (2)

Under these set of coordinates the relevant ($r-t$)-sector of the metric (1)
takes the form,

$\displaystyle ds^{2}=\frac{F(r)}{2}(du~{}dv+dv~{}du)$ (3)

Chiral conditions, to be discussed in the next section, are most appropriately
described in these coordinates.

Chirality conditions: Consider the Klein-Gordon equation for a massless scalar
particle governed by the metric (3),

$\displaystyle
g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi=\frac{4}{F}\partial_{u}\partial_{v}\phi=0.$
(4)

The general solution of this can be taken as
$\phi(u,v)=\phi^{(R)}(u)+\phi^{(L)}(v)$ where $\phi^{(R)}(u)$ and
$\phi^{(L)}(v)$ are the right (outgoing) and left (ingoing) modes satisfying

$\displaystyle\nabla_{v}\phi^{(R)}=0,\,\,\nabla_{u}\phi^{(R)}\neq 0;\,\,\,\
\nabla_{u}\phi^{(L)}=0,\,\,\,\nabla_{v}\phi^{(L)}\neq 0.$ (5)

These equations are expressed simultaneously as,

$\displaystyle\nabla_{\mu}\phi=\pm\bar{\epsilon}_{\mu\nu}\nabla^{\nu}\phi=\pm\sqrt{-g}\epsilon_{\mu\nu}\nabla^{\nu}\phi$
(6)

where $+(-)$ stand for left (right) mode and $\epsilon_{\mu\nu}$ is the
numerical antisymmetric tensor with $\epsilon_{uv}=\epsilon_{tr}=-1$. This is
the chirality condition 111In analogy with studies in 2d CFT this condition is
usually referred as holomorphy condition.. In fact the condition (6) holds for
any chiral vector $J_{\mu}$ in which case
$J_{\mu}=\pm{\bar{\epsilon}}_{\mu\nu}J^{\nu}$. Likewise, the chirality
condition for the energy-momentum tensor is [17],

$\displaystyle
T_{\mu\nu}=\pm\frac{1}{2}(\bar{\epsilon}_{\mu\sigma}T^{\sigma}_{\nu}+\bar{\epsilon}_{\nu\sigma}T^{\sigma}_{\mu})+\frac{1}{2}g_{\mu\nu}T^{\alpha}_{\alpha}$
(7)

This shows that corresponding to right mode (i.e. ($-$) sign)

$\displaystyle T^{(R)}_{vv}=0,\,\,\,\,T^{(R)}_{uu}\neq 0$ (8)

which is analogous to the first equation of (5). The other sign of (7) gives
the energy-momentum tensor for the left mode. This essentially corresponds to
the interchange $u\leftrightarrow v$ and $L\leftrightarrow R$. In the next
section, using these chirality conditions we will derive the explicit form for
the gravitational anomaly that reproduces the Hawking flux.

Chirality, gravitational anomaly and Hawking flux: It is well known that for a
non-chiral (vector like) theory it is not possible to simultaneously preserve,
at the quantum level, general coordinate invariance as well as conformal
invariance. Since the former invariance is fundamental in general relativity,
conformal invariance is sacrificed leading to a nonvanishing trace of the
stress tensor, called the trace anomaly. Using this trace anomaly and the
chirality condition we will derive an expression for the chiral gravitational
(diffeomophism) anomaly from which the Hawking flux is computed.

The energy-momentum tensor near an evaporating black hole is split into a
traceful and traceless part by [18],

$\displaystyle T_{\mu\nu}=\frac{R}{48\pi}g_{\mu\nu}+\theta_{\mu\nu}$ (9)

where $\theta_{\mu\nu}$ is symmetric (i.e. $\theta_{\mu\nu}=\theta_{\nu\mu}$),
so that it preserves the symmetricity of $T_{\mu\nu}$, and traceless (i.e.
$\theta_{\mu}^{\mu}=0$ so that in $u,v$ coordinates $\theta_{uv}=0$). The
traceful part is contained in the first piece leading to the trace anomaly,
$T^{\mu}_{\mu}=\frac{R}{24\pi}$. Also, since general coordinate invariance is
preserved, $\nabla^{\mu}T_{\mu\nu}=0$, from which it follows that the
solutions of $\theta_{\mu\nu}$ satisfy,

$\displaystyle\nabla^{\mu}\theta_{\mu\nu}=-\frac{1}{48\pi}\nabla_{\nu}R$ (10)

Now the energy-momentum tensor (9) can be regarded as the sum of the
contributions from the right and left moving modes. Symmetry principle tells
that the contribution from one mode is exactly equal to that from the other
mode, only that $u,v$ have to be interchanged. Since $T_{\mu\nu}$ is symmetric
we have $T_{\mu\nu}=T_{\mu\nu}^{(R)}+T_{\mu\nu}^{(L)}$ with

$\displaystyle
T_{\mu\nu}^{(R/L)}=\frac{R}{96\pi}g_{\mu\nu}+\theta_{\mu\nu}^{(R/L)}$ (11)

where $\theta_{\mu\nu}=\theta_{\mu\nu}^{(R)}+\theta_{\mu\nu}^{(L)}$ (in
analogy with $T_{\mu\nu}$). Therefore the chirality condition (8) and the
traceless condition of $\theta_{\mu\nu}$ immediately show

$\displaystyle\theta_{uv}^{(R)}=0,\,\,\,\theta_{vv}^{(R)}=0,\,\,\,\,\theta_{uu}^{(R)}\neq
0;\,\,\,\theta_{uv}^{(L)}=0,\,\,\,\theta_{uu}^{(L)}=0,\,\,\,\,\theta_{vv}^{(L)}\neq
0$ (12)

The trace anomaly for the chiral components follows from (11) and (12),

$\displaystyle
T{{}^{\mu}_{\mu}}{{}^{(R)}}=T{{}^{\mu}_{\mu}}{{}^{(L)}}=\frac{1}{2}T^{\mu}_{\mu}=\frac{R}{48\pi}$
(13)

To find out the diffeomorphism anomaly for the chiral components we will use
(11). Considering only the right mode, for example, we have

$\displaystyle\nabla^{\mu}T_{\mu\nu}^{(R)}=\frac{1}{96\pi}\nabla_{\nu}R+\nabla^{\mu}\theta_{\mu\nu}^{(R)}$
(14)

Next, using (10) and (12) for the right mode we obtain,

$\displaystyle\nabla^{\mu}\theta_{\mu
u}^{(R)}=-\frac{1}{48\pi}\nabla_{u}R;\,\,\,\nabla^{\mu}\theta_{\mu v}^{(R)}=0$
(15)

Substituting these in (14) we get, once for $\nu=u$ and then $\nu=v$,

$\displaystyle\nabla^{\mu}T_{\mu
u}^{(R)}=-\frac{1}{96\pi}\nabla_{u}R;\,\,\,\,\,\nabla^{\mu}T_{\mu
v}^{(R)}=\frac{1}{96\pi}\nabla_{v}R.$ (16)

Therefore, combining both the above results yields

$\displaystyle\nabla^{\mu}T_{\mu\nu}^{(R)}=\frac{1}{96\pi}\bar{\epsilon}_{\nu\lambda}\nabla^{\lambda}R$
(17)

which is the chiral (gravitational) anomaly for the right mode. Similarly the
chiral anomaly for left mode can also be obtained which has a similar form
except for a minus sign on the right side of (17). This anomaly is in
covariant form and so it is also called the covariant gravitational anomaly.
The structure, including the normalization, agrees with that found by using
explicit regularization of the chiral stress tensor [19, 20].

From (17) and (13) a simple relation follows between the gravitational anomaly
(${\cal{A}}_{\nu}$) and the trace anomaly ($T$),

$\displaystyle{\cal{A}}_{\nu}=\frac{1}{2}\bar{\epsilon}_{\nu\lambda}\nabla^{\lambda}T.$
(18)

Such a relation is not totally unexpected since covariant expressions must
involve the Ricci scalar. However (18) should not be interpreted as a Wess-
Zumino consistency condition which involves only ‘consistent’ expressions
[20]. Here, on the contrary, we are dealing with covariant expressions.

The covariant anomaly (17) is now used to obtain the Hawking flux. As was
mentioned earlier the effective two dimensional theory near the horizon
becomes chiral. The chiral theory has the anomaly (17). Taking its $\nu=u$
component we obtain,

$\displaystyle\partial_{r}T_{uu}^{(R)}=\frac{F}{96\pi}\partial_{r}R=\frac{F}{96\pi}\partial_{r}(F^{\prime\prime})=\frac{1}{96\pi}\partial_{r}(FF^{\prime\prime}-\frac{F^{\prime
2}}{2})$ (19)

which yields,

$\displaystyle
T_{uu}^{(R)}=\frac{1}{96\pi}\Big{(}FF^{{}^{\prime\prime}}-\frac{F^{\prime
2}}{2}\Big{)}+C_{uu}$ (20)

where $C_{uu}$ is an integration constant. Imposing the usual boundary
condition which requires $T_{uu}^{(R)}(r\rightarrow r_{H})=0$, implying that a
freely falling observer sees a finite amount of flux at the outer horizon,
leads to $C_{uu}=\frac{F^{\prime 2}(r_{H})}{192\pi}$. This condition on the
outgoing modes is similar to that of the Unruh vacuum [17]. The corresponding
condition on the ingoing mode for the Unruh vacuum is satisfied by default
since, due to chirality, these are absent ($T_{vv}^{(R)}=0$). Note, however,
that the Unruh condition on the ingoing modes
$T_{vv}^{(R)}(r\rightarrow\infty)=0$ is applied at asymptotic infinity where
the theory is non-chiral. This does not affect our interpretation since,
asymptotically, the anomaly (17) vanishes. Hence the results from the chiral
expressions will agree with the non-chiral ones at asymptotic infinity.
Indeed, the Hawking flux, obtained by taking the asymptotic infinity limit
($r\rightarrow\infty$) of (20) [6, 7, 8, 11, 13],

$\displaystyle T_{uu}^{(R)}(r\rightarrow\infty)=C_{uu}=\frac{F^{\prime
2}(r_{H})}{192\pi}=\frac{K^{2}}{48\pi}$ (21)

where $K=\frac{F^{\prime}(r_{H})}{2}$ is the surface gravity of the black
hole, reproduces the known result corresponding to the Hawking temperature
$T_{H}=\frac{K}{2\pi}$. The other terms in (20) drop out due to asymptotic
flatness.

Chirality, quantum tunneling and Hawking temperature: Here, using the
chirality condition (6), we will derive the tunneling probability, which will
eventually yield the Hawking temperature. Under the metric (1) this condition
corresponds to,

$\displaystyle\partial_{t}\phi(r,t)=\pm F(r)\partial_{r}\phi(r,t)$ (22)

As before $+(-)$ stand for left (right) mode. Putting the standard WKB ansatz

$\displaystyle\phi(r,t)=e^{-\frac{i}{\hbar}S(r,t)}$ (23)

in (22), where $S(r,t)$ is the action, we get the familiar semiclassical
Hamilton-Jacobi equation

$\displaystyle\partial_{t}S(r,t)=\pm F(r)\partial_{r}S(r,t)$ (24)

which is the basic equation in the tunneling mechanism for studying Hawking
radiation. This has been derived earlier from the Klein-Gordon equation with
the background metric (1) and the ansatz (23) [4, 10].

Now since the metric (1) is stationary we choose a solution for $S(r,t)$ as

$\displaystyle S(r,t)=\omega t+S(r)$ (25)

where $\omega$ is the energy of the particle. Substituting this in (24) a
solution for $S(r)$ is obtained. Inserting this back in (25) yields,

$\displaystyle S(r,t)=\omega t\pm\omega\int\frac{dr}{F(r)}$ (26)

It is important to note that expressing (26) in the null tortoise coordinates
(see (2)) defined inside and outside of the event horizon and then
substituting in (23) one can obtain the right and left modes for both sectors:

$\displaystyle\Big{(}\phi^{(R)}\Big{)}_{\textrm{in}}=e^{-\frac{i}{\hbar}\omega
u_{\textrm{in}}};\,\,\,\Big{(}\phi^{(L)}\Big{)}_{\textrm{in}}=e^{-\frac{i}{\hbar}\omega
v_{\textrm{in}}}$
$\displaystyle\Big{(}\phi^{(R)}\Big{)}_{\textrm{out}}=e^{-\frac{i}{\hbar}\omega
u_{\textrm{out}}};\,\,\,\Big{(}\phi^{(L)}\Big{)}_{\textrm{out}}=e^{-\frac{i}{\hbar}\omega
v_{\textrm{out}}}$ (27)

which satisfy the condition (5). Precisely these modes were used previously to
find the trace anomaly [18] as well as the chiral (gravitational) anomaly [19]
by the point splitting regularization technique. In our formulation these
modes (27) are a natural consequence of chirality.

Now in the tunneling formalism a virtual pair of particles is produced in the
black hole. One of this pair can quantum mechanically tunnel through the
horizon. This particle is observed at infinity while the other goes towards
the center of the black hole. While crossing the horizon the nature of the
coordinates changes. This can be explained in the following way. The Kruskal
time ($T$) and space ($X$) coordinates inside and outside the horizon are
defined as,

$\displaystyle
T_{\textrm{in}}=e^{K(r_{*})_{\textrm{in}}}~{}{\textrm{cosh}}(Kt_{\textrm{in}});\,\,\,X_{\textrm{in}}=e^{K(r_{*})_{\textrm{in}}}~{}{\textrm{sinh}}(Kt_{\textrm{in}})$
$\displaystyle
T_{\textrm{out}}=e^{K(r_{*})_{\textrm{out}}}~{}{\textrm{sinh}}(Kt_{\textrm{out}});\,\,\,X_{\textrm{out}}=e^{K(r_{*})_{\textrm{out}}}~{}{\textrm{cosh}}(Kt_{\textrm{out}})$
(28)

where, as before, $K=\frac{F^{\prime}(r_{H})}{2}$ is the surface gravity of
the black hole. These two sets of coordinates are connected by the relations,

$\displaystyle
t_{\textrm{in}}=t_{\textrm{out}}-i\frac{\pi}{2K};\,\,\,\,(r_{*})_{\textrm{in}}=(r_{*})_{\textrm{out}}+i\frac{\pi}{2K}$
(29)

Therefore following the definition (2) we obtain the relations connecting the
null coordinates defined inside and outside the horizon,

$\displaystyle
u_{\textrm{in}}=t_{\textrm{in}}-(r_{*})_{\textrm{in}}=u_{\textrm{out}}-i\frac{\pi}{K}$
$\displaystyle
v_{\textrm{in}}=t_{\textrm{in}}+(r_{*})_{\textrm{in}}=v_{\textrm{out}}$ (30)

Since the left moving mode travels towards the center of the black hole, its
probability to go inside is

$\displaystyle P^{(L)}=|\phi^{(L)}_{\textrm{in}}|^{2}=1$ (31)

i.e. the left moving (ingoing) mode is trapped inside the black hole. On the
other hand the right moving mode ($\phi^{(R)}_{\textrm{in}}$) tunnels through
the event horizon. To calculate the tunneling probability as seen by an
external observer it is first necessary to rewrite $u_{\textrm{in}}$ in terms
of $u_{\textrm{out}}$ using(30). We find, from (27),

$\displaystyle
P^{(R)}=|\phi^{(R)}_{\textrm{in}}|^{2}=|e^{-\frac{i}{\hbar}\omega(u_{\textrm{out}}-i\frac{\pi}{K})}|^{2}=e^{-\frac{2\pi\omega}{\hbar
K}}$ (32)

Then using the principle of “detailed balance” [4, 10]
$P^{(R)}=e^{-\frac{\omega}{T_{H}}}P^{(L)}=e^{-\frac{\omega}{T_{H}}}$ yields
the Hawking temperature (in the units of $\hbar=1$) as,

$\displaystyle T_{H}=\frac{K}{2\pi}$ (33)

This is the standard Hawking temperature corresponding to the flux (21).

As we observe the ingoing modes are trapped and do not play any role in the
computation of the Hawking temperature. A similar feature occurs in the
anomaly approach where the ingoing modes are neglected leading to a chiral
theory that eventually yields the flux. These observations provide a physical
picture of chirality connecting the tunneling and anomaly methods.

Conclusions: We have shown that the notion of chirality pervades the anomaly
and tunneling formalisms thereby providing a close connection between them.
This is true both from a physical as well as algebraic perspective. The chiral
restrictions play a pivotal role in the abstraction of the anomaly from which
the flux is computed. The same restrictions, in the tunneling formalism, lead
to the Hawking temperature corresponding to that flux.

A dimensional reduction is known to reduce the theory effectively to a two
dimensional conformal theory near the event horizon. The ingoing (left moving)
modes are lost inside the horizon. They cannot contribute to the near horizon
theory thereby rendering it chiral and, hence, anomalous. Using the
restrictions imposed by chirality we obtained a form for this (gravitational)
anomaly, manifested by a nonconservation of the stress tensor, by starting
from the familiar form of the trace anomaly. From a knowledge of the
gravitational anomaly we were able to obtain the flux.

The chirality constraints were then exploited to obtain the equations for the
ingoing and outgoing modes in the tunneling formalism, following the standard
geometrical (WKB) approximation. We reformulated the tunneling mechanism to
highlight the role of coordinate systems in the chiral framework. A specific
feature of this reformulation is that explicit treatment of the singularity in
(26) is not required since we do not carry out the integration. Only the modes
inside ($\phi_{\textrm{in}}$) and outside ($\phi_{\textrm{out}}$) the horizon,
both of which are well defined, are required. The singularity now manifests in
the complex transformations (29) that connects these modes across the horizon.
The probability for finding the ingoing modes was shown to be unity. These
modes do not play any role in the tunneling approach which is the exact
analogue of omitting them when considering the effective near horizon theory
in the anomaly method.

It is useful to observe that the crucial role of chirality in both approaches
is manifested in the near horizon regime. This reaffirms the universality of
the Hawking effect being governed by the properties of the event horizon.

## References

  * [1] S.W.Hawking, Commun. Math. Phys. 43, 199 (1975).
  * [2] G.W.Gibbons and S.W.Hawking, Phys. Rev. D 15, 2752 (1977).
  * [3] S.M.Christensen and S.A.Fulling, Phys. Rev. D 15, 2088 (1977).
  * [4] K.Srinivasan and T.Padmanabhan, Phys. Rev. D 60, 024007 (1999) [arXiv:gr-qc/9812028].
  * [5] M.K.Parikh and F.Wilczek, Phys. Rev. Lett. 85, 5042 (2000) [arXiv:hep-th/9907001].
  * [6] S.P.Robinson and F.Wilczek, Phys. Rev. Lett. 95, 011303 (2005) [arXiv:gr-qc/0502074].
S.Iso, H.Umestsu and F.Wilczek, Phys. Rev. Lett. 96, 151302 (2006) [arXiv:hep-
th/0602146].

  * [7] R.Banerjee and S.Kulkarni, Phys. Rev. D 77, 024018 (2008) [arXiv:0707.2449].
  * [8] R.Banerjee and S.Kulkarni, Phys. Lett. B 659, 827 (2008) [arXiv:0709.3916].
  * [9] M.Arzano, A.J.M.Medved and E.C.Vagenas, JHEP 0509, 037 (2005) [arXiv:hep-th/0505266].
E.T.Akhmedov, V.A.Akhmedova and D.Singleton, Phys. Lett. B 642, 124 (2006)
[arXiv:hep-th/0608098].
M.Angheben, M.Nadalini, L.Vanzo and S.Zerbini, JHEP 0505, 014 (2005)
[arXiv:hep-th/0503081].
R.Kerner and R.B.Mann, Phys. Rev. D 73, 104010 (2006) [arXiv:gr-qc/0603019].
P.Mitra, Phys. Lett. B 648, 240 (2007) [arXiv:hep-th/0611265].
R.Banerjee and B.R.Majhi, Phys. Lett. B 662, 62 (2008) [arXiv:0801.0200].

  * [10] R.Banerjee and B.R.Majhi, JHEP 0806, 095 (2008) [arXiv:0805.2220].
  * [11] R.Banerjee, [arXiv:0807.4637].
  * [12] S.Ghosh, [arXiv:0804.2999].
  * [13] S.Iso, T.Morita and H.Umestsu, Phys. Rev. D 77, 045007 (2008) [arXiv:0710.0456].
L.Bonora and M.Cvitan, JHEP 0805, 071 (2008)[arXiv:0804.0198].
L.Bonora, M.Cvitan, S.Pallua and I.Smolic, [arXiv:0808.2360].

  * [14] T.Morita, [arXiv:0811.1741].
  * [15] T.Damour and R.Ruffini, Phys. Rev. D 14, 332 (1976).
  * [16] S.Carlip, Phys. Rev. Lett. 82, 2828 (1999) [arXiv:hep-th/9812013].
  * [17] R.Banerjee and S.Kulkarni, [arXiv:0810.5683].
  * [18] P.C.W.Davies and S.A.Fulling, Proc. R. Soc. Lond. A 354, 59 (1977).
  * [19] S.A.Fulling, General Relativity and Gravitation 18, 609 (1986).
  * [20] R.Bertlmann and E.Kohlprath, Ann. Phys. 288, 137 (2001).


$\underset{}{}$Sequentially estimating the dynamic contact angle
of$\underset{}{}$
---
$\underset{}{}$sessile saliva droplets in view of SARS-CoV-2$\underset{}{}$

Sudeep R. Bapat
---
$\underset{}{}$Indian Institute of Management, Indore, India$\underset{}{}$
E-mail<EMAIL_ADDRESS>

Abstract: Estimating the contact angle of a virus infected saliva droplet is
seen to be an important area of research as it presents an idea about the
drying time of the respective droplet and in turn of the growth of the
underlying pandemic. In this paper we extend the data presented by Balusamy,
Banerjee and Sahu [“Lifetime of sessile saliva droplets in the context of
SARS-CoV-2,” Int. J. Heat Mass Transf. 123, 105178 (2021)], where the contact
angles are fitted using a newly proposed half-circular wrapped-exponential
model, and a sequential confidence interval estimation approach is established
which largely reduces both time and cost with regards to data collection.

Keywords: SARS-CoV-2; COVID-19; Circular data; Wrapped exponential; Contact
angle

1\. INTRODUCTION

SARS-CoV-2 (virus which causes COVID-19) has severely impacted more than 200
countries worldwide, with over 180 million cases until the end of June, 2021.
The span of this virus was so fast and devastating that the World Health
Organization declared the outbreak as a Public Health Emergency of
International Concern on 30 January, 2020, whereas a global pandemic on 11
March, 2020. Spreading of such respiratory diseases is largely caused due to
respiratory droplets of saliva (of an infected person) during coughing,
sneezing or even moist speaking. A recent reference paper in this regard is by
Balusamy et al. [1]. Understanding the lifetime of such droplets is hence an
important area of research, which could be handled by studying the fluid
dynamics of such droplets in air. One may refer to Mittal et al. [2] who
analyze the flow-physics of virus laden respiratory droplets, or Bhardwaj and
Aggarwal [3] who analyze the likelihood of survival of a virus laden droplet
on a solid surface. Further, it has been studied that such respiratory
droplets have a tendency to increase their lifetime on coming in contact with
a surface based on its properties. Vejerano and Marr [4] studied the physico-
chemical characteristics of evaporating respiratory fluid droplets and found
out that a typical saliva droplet also contains NaCl, mucin (protein) and a
certain surfactant in fixed amounts. In addition to the droplet composition,
the evaporation rate of a droplet also depends on environmental conditions and
factors such as temperature, relative humidity, droplet volume and the contact
angle which the droplet makes with the surface. A specific analysis was
carried out in [3] where the authors examined the drying time of a deposited
droplet in two different temperatures namely, $25^{\circ}$C and $40^{\circ}$C
which represent an air-conditioned room and a summer afternoon respectively.
The contact angle and humidity were set at $30^{\circ}$ and $50\%$. Studying
the drying time of a droplet plays an important role as it well related to the
the survival of the droplet and in turn to the growth of the pandemic.
Chaudhuri et al. [5] tested this hypothesis using suspended droplets in air,
whereas [3] compared the growth of infection with the drying time in different
cities globally. They verified that for a $5$ nL droplet, a higher drying time
corresponds to a higher pandemic growth rate. Hence, when a droplet evaporates
slowly, the chance of the survival of the virus is enhanced.

Specifically, the initial contact angle, which measures the angle that a
droplet makes with the surface plays a big role in determining the lifetime of
it. Different contact angles are predominant with different surfaces i.e.,
droplets on glass, wood, stainless steel, cotton or the touchscreen of a
smartphone tend to make angles varying from $5^{\circ}$ to $95^{\circ}$. It is
also intuitive that a contact angle cannot exceed $180^{\circ}$. Figure 1
contains pictorial representations of two different droplets making different
angles with the surface. The left image shows a water droplet on cloth, making
a high contact angle due to the hydrophobic property of the cloth. Whereas the
image on the right shows a water droplet on a lotus leaf, again making a high
contact angle. Both the images are borrowed from Wikipedia under the license
CC BY-SA 3.0.

|
---|---
(a) water drop on cloth | (b) water drop on a lotus leaf
Figure 1. Water droplets making contact angles greater than $90^{\circ}$ on
two different surfaces. Both the
images are borrowed from Wikipedia under the license CC BY-SA 3.0

A dynamic contact angle is the one which is measured as the droplet changes
its size as it moves quickly over the surface. One may again refer to [3] or
[1] for more details. However it is also true that measuring such contact
angles (initial or dynamic) involves a lot of struggle and cost, as it has to
be carried out using heavy apparatus. Some of the existing methods for contact
angle determination include the sessile droplet method, where the angle is
measured using a “contact angle goniometer”, the pendant drop method which is
used to measure angles for pendant drops, the dynamic sessile drop method
which is similar to a sessile drop method but requires the drop to be modified
or a single-fiber meniscus method where the shape of the meniscus on the fiber
is directly imaged through a high resolution camera. One may refer to Albert
et al. [6] for an overview of other techniques. Hence to estimate a dynamic
contact angle of a droplet, a reduction in the number of observations required
to carry out the estimation is highly beneficial. In this paper, we thus
introduce a sequential estimation technique which is outlined in Section 2.

Now since the aim of this paper is to estimate a certain contact angle, it
makes more sense to apply a circular model rather than a usual linear one on
the concerned data. Literature on such models is vast and ever expanding. A
few other examples where a circular model is appropriate involve orientations
of the major axis of termite mounds, the angles of slope of different
sedimentary layers of an exposed rock-face or the walking directions of long
legged desert ants etc. In all these examples, the observations are either
certain directions, or angles measured in degrees or radians. Such
observations are often measured either clockwise or counter-clockwise from
some reference direction, usually called as the zero direction. Over years, a
usual technique to design new circular distributions is to wrap a linear
distribution over a full circle. However as seen before, since the contact
angles of any droplet is necessarily less than $180^{\circ}$, an adjusted
model which is capable of taking values only on half-a-circle seems more
appropriate. In this context, we introduce a new model called as the half-
circular wrapped-exponential distribution to model our data. In general, a few
notable books covering circular models which one can refer to are by Mardia
and Jupp [7], Rao and Sengupta [8], Rao and Girija [9] or Ley and Verdebout
[10], among others.

2\. Data modeling and analysis

The particular dataset analyzed for this experiment is a pseudo dataset which
is an extended version of the one borrowed from [1] and consists of the
temporal variations of the dynamic contact angles in degrees (simply called as
contact angles from now on) of the droplet normalized with the initial contact
angle, $\theta/\theta_{0}$. The particular setting used for this experiment is
as follows: the relative humidity (RH) is controlled at $50\%$, the initial
droplet volume $(V_{0})$ is $10$ nL, the molality of the saliva $(M)$ is
$0.154$ mol/kg, temperature $(T)$ is $30^{\circ}$, the surfactant parameter
$(\Psi)$ is $10$ and the initial contact angle $(\theta_{0})$ is $50^{\circ}$.
One may refer to Figure 2a in [1] for a pictorial description of the dataset.
As there was not an access to the actual observations, we adopted the
following approach: for brevity alone, we only focused on the curve
representing RH$=50\%$. Using an online tool, we extracted the $(x,y)$
coordinates for each of its 20 observations. We converted these normalized
contact angles to actual contact angles $(\theta)$ and finally translated
those into radians. Table 1 lists all these observations for convenience. Now,
to include more observations in the analysis, we first assumed a linear
relationship between “time” and “contact angles” (CA), fitted several
polynomial regression models and picked the following third-order model which
fitted better with a $R^{2}$ value of $0.9613$.

$CA=0.985-8.45\times 10^{-3}\,time+2.34\times 10^{-5}\,time^{2}-2.05\times
10^{-8}\,time^{3}$ (1)

Figure 2 contains a scatterplot of the raw data (1a) and the fitted polynomial
regression model superimposed on it (1b).

Table 1. Extracted dataset containing the temporal variations
---
of the contact angles (in radians)
$\underset{}{\overset{}{\text{Time (sec)}}}$ | CA | Time (sec) | CA | Time (sec) | CA
$10$ | $0.811$ | $88.75$ | $0.379$ | $287.5$ | $0.034$
$25$ | $0.794$ | $100$ | $0.261$ | $325$ | $0.031$
$55$ | $0.689$ | $118.75$ | $0.218$ | $381.25$ | $0.028$
$58.75$ | $0.654$ | $137.5$ | $0.157$ | $437.5$ | $0.026$
$66.25$ | $0.593$ | $175$ | $0.109$ | $493.75$ | $0.023$
$77.25$ | $0.471$ | $212.5$ | $0.052$ | $550$ | $0.020$
$83.15$ | $0.436$ | $250$ | $0.035$ |  |

|
---|---
(a) plot of the raw data | (b) superimposed polynomial model
Figure 2. Temporal variations of the contact angles

We then assumed a vector of times ranging over $5-300$ seconds with a jump of
$1$ second in between, and predicted the contact angles according to the above
model. Thus, our final pseudo dataset consists of $296$ observations according
to our construction.

2.1. A half-circular wrapped-exponential model for the contact angles

For a start, Figure 3 shows a pictorial distribution of our pseudo data placed
around a circle. Purposefully, we have stacked the closely lying observations
for a better visualization and as one can observe, all the observations lie
entirely between $0$ and $\pi/2$ radians. As seen before, wrapping a linear
density over a circle is a suitable choice to model such observations. In this
case, since the linear curve seen in Figure 2 shows an exponential decline, it
makes sense to choose some of the lifetime distributions and wrapping them
around a circle. Now as discussed before, since any contact angle of a droplet
is always less than $\pi$ radians it makes more sense to fit a distribution
which takes values only on a semicircle. In literature, not many such
distributions have been proposed. One such example is of a half-circular
distribution which was introduced by Rambli et al. [11], who converted a Gamma
distribution to a half-circular one and fitted it to the angle which measures
the posterior corneal curvature of an eye. In a similar spirit, we now
introduce a half-circular wrapped-exponential (HCWE) distribution with
parameter $\lambda$. An intuitive construction is through the following
transformation: $X_{w}=X(\text{mod}\,\pi)$, where $X$ is a linear exponential
random variable with pdf $f(x)=\lambda e^{-\lambda x},\,x>0,\lambda>0$.
Interestingly, another easy construction is to simply truncate $X$ over
$[0,\pi)$. Its pdf, cdf and characteristic functions are as follows,

$f_{w}(\theta)=\frac{\lambda e^{-\lambda\theta}}{1-e^{-\pi\lambda}},\hskip
8.5359pt\theta\in[0,\pi)$ (2)

$F_{w}(\theta)=\frac{1-e^{-\lambda\theta}}{1-e^{-\pi\lambda}},\hskip
8.5359pt\theta\in[0,\pi)$ (3)

$\phi_{p}=\frac{1}{1-ip/\lambda},\hskip 8.5359ptp=0,\pm 1,\pm 2,...$ (4)

Consequently, the mean direction happens to be,

$\mu_{0}=\tan^{-1}\frac{1}{\lambda},\hskip 8.5359pt\lambda>0$ (5)

Now for a comparison, we tried to fit several other wrapped distributions to
the data namely, the wrapped-exponential by Jammalamadaka and Kozubowski [12],
transmuted wrapped-exponential by Yilmaz and Biçer [13] and wrapped-Lindley by
Joshi and Jose [14]. For completeness, we also fit a von-Mises distribution
which is one of the widely used circular models. Table 2 contains the log-
likelihood values and the AICs for these five models. As one can observe, the
half-circular wrapped-exponential model fits better than the others. It is
also seen to be a significant fit with a p-value of $0.18$ using the
Kolmogorov-Smirnov test, and the estimated $\lambda$ value equals $3.69$. On
using Eq. (5), the estimated mean direction equals $0.2646$ radians. Figure 4
contains a set of goodness of fit plots for the HCWE$(\lambda)$ distribution.
All these fits and plots were carried out using the $``circular"$ and
$``fitdistrplus"$ packages in $R$.

Figure 3. Raw circular plot of the pseudo data

Figure 4. Goodness of fit plots for the HCWE$(\lambda)$ density on the pseudo
data

Now, since in practice the value of $\lambda$ will be unknown, we develop a
sequential fixed-width confidence interval to estimate $\lambda$ which in turn
will give us an estimate for the mean direction $\mu_{0}$ of the contact
angle, which will give us a fair idea about the drying time of the droplet.

Table 2. Comparing model fits to the pseudo data
---
$\underset{}{\overset{}{\text{Model}}}$ | Log-likelihood | AIC
von-mises | $-36.72$ | $77.44$
wrapped-exponential | $90.75$ | $-179.50$
half-circular wrapped-exponential | 92.56 | -181.52
transmuted wrapped-exponential | $-7.01$ | $18.02$
wrapped-Lindley | $89.94$ | $-177.88$

2.2. A sequential fixed-width confidence interval

In general, a sequential rule consists of identifying a stopping variable,
which determines the optimal sample size to be used in the experiment. This
technique largely reduces the number of observations needed for the inference
part, which proves to be beneficial as it reduces both time and cost.
Literature on sequential estimation methodologies is vast and still being
explored. In particular, a few recent works aimed at finding appropriate
confidence intervals include, Banerjee and Mukhopadhyay [15], who developed a
general sequential fixed-accuracy confidence interval, Mukhopadhyay and
Banerjee [16], who looked at constructing bounded length intervals, Bapat [17,
18] who constructed fixed-accuracy intervals for parameters under an inverse
Gaussian and bivariate exponential models or Khalifeh et al. [19], who derived
fixed-accuracy intervals for the reliability parameter of an exponential
distribution.

To summarize, a fixed-width interval (FWI) aims at simultaneously controlling
the width of the interval (say, $d)$ and the confidence limit $(1-\alpha)$.
Such an interval is clearly symmetric around the parameter. It turns out that
there does not exist any fixed sample size procedure to tackle this problem
and one has to resort to a sequential setup. However, a certain drawback of
this method is: even though a parameter is entirely positive, the lower bound
of a FWI can assume negative values. A fix to this is to construct a fixed-
accuracy interval (FAI), which assumes a fixed-accuracy value (say, $d)$. A
FAI happens to be symmetric around $\log$ of the parameter. An introductory
paper to this approach is Banerjee and Mukhopadhyay [20]. Even in this case it
may happen, that if the parameter space is bounded (say from above by $U$), a
FAI may contain bounds which cross $U$. Hence, [15] came up with a bounded-
length fixed-accuracy interval (BLFAI) as a fix. In our case, we aim at
constructing a fixed-width interval as outlined next.

Let $\theta_{1},\theta_{2},...$ be the dynamic contact angles of a droplet,
measured using a suitable technique. Then, for some pre-fixed width $d$, a
confidence interval of $\lambda$ takes the following form,

$I_{n}=\left\\{\lambda:\lambda\in\left[\hat{\lambda}_{n}-d,\hat{\lambda}_{n}+d\right]\right\\},$
(6)

where $\hat{\lambda}_{n}$ is the MLE of $\lambda$, which is consistent and
asymptotically normal with the following representation,

$\sqrt{n}(\hat{\lambda}_{n}-\lambda)\xrightarrow{D}N\left(0,\sigma^{2}_{\hat{\lambda}_{n}}\right),$
(7)

where $\sigma^{2}_{\hat{\lambda}_{n}}$ is the variance of the MLE and
$\xrightarrow{D}$ stands for convergence in distribution. Now, for $I_{n}$ to
include $\lambda$ with a pre-fixed coverage probability $1-\alpha$, the
required fixed sample size can be found out as follows,

$P\left(\hat{\lambda}_{n}-d\leq\lambda\leq\hat{\lambda}_{n}+d\right)=1-\alpha$
$\Rightarrow n^{*}\equiv
n_{d}^{*}=\left(\frac{z_{\alpha/2}}{d}\right)^{2}\sigma^{2}_{\hat{\lambda}_{n}},$
(8)

where $z_{\alpha/2}$ is the upper $100(\alpha/2)\%$ point of a standard normal
distribution. Since $n^{*}$ is an unknown quantity, we now propose the
following sequential methodology:

We first fix an integer $m(>1)$ often called as the “pilot sample size” and
obtain a pilot sample $\theta_{1},\theta_{2},...,\theta_{m}$ from a
HCWE$(\lambda)$ density as given in (5). We then aim to collect an additional
observation at every stage, until sampling is terminated according to the
following stopping rule:

$N=\inf\left\\{n\geq
m:n\geq\left(\frac{z_{\alpha/2}}{d}\right)^{2}\hat{\sigma}^{2}_{\hat{\lambda}_{n}}\right\\},$
(9)

where $\hat{\sigma}^{2}_{\hat{\lambda}_{n}}$ is the estimated variance of the
MLE. We then have a final set of observations
$\theta_{1},\theta_{2},...,\theta_{N}$ and will estimate $\lambda$ using the
interval,

$I_{N}=\left[\hat{\lambda}_{N}-d,\hat{\lambda}_{N}+d\right]=\left[L_{N},U_{N}\right]\,\text{(say)}.$
(10)

The stopping variable $N$ from (12) follows properties such as asymptotic
first-order efficiency and asymptotic consistency. We leave out the proofs for
brevity. One may refer to Theorem 3.1 of [17]. Finally, we estimate the mean
direction $\mu_{0}$ using an interval,

$J_{N}=\left[\tan^{-1}\frac{1}{U_{N}},\,\,\tan^{-1}\frac{1}{L_{N}}\right].$
(11)

We now outline a stepwise procedure to tackle a practical problem through the
above methodology.

  * Step 1:

For a certain specific liquid droplet, observe the contact angles over equally
spaced time intervals and note down the first $m$ angles
$(\theta_{1},\theta_{2},...,\theta_{m})$ over the first $m$ time points
$t_{1},t_{2},...,t_{m}.$

  * Step 2:

After $t_{m}$, collect observations (i.e. observe contact angles) one-at-a-
time according to the stopping rule given in (9).

  * Step 3:

Once the stopping rule is executed, observe the value of $N$, find out an
interval for $\lambda$ as per (10) and ultimately find a subsequent interval
for the mean direction $\mu_{0}$ according to (11).

  * Step 4:

Using the interval for $\mu_{0}$, find out a rough interval for the average
drying time of the droplet by predicting using the following inverted
polynomial regression model $(R^{2}=0.98)$ (i.e. by assuming “time” as the
response and “contact angle” as the predictor.

$time=266.96-872.293\,CA+1329.892\,CA^{2}-763.05\,CA^{3}$ (12)

Hence, for our complete pseudo data,
$\hat{\lambda}_{n}=3.69,\hat{\mu}_{0}=0.2646$ and the estimated drying time
equals $115.13$ seconds. We now apply the above procedure to our observed
pseudo data with a small adjustment: we first randomize the entire data,
sample $250$ observations, and sort them. This kind of an approach gives a
good representation of the actual data in every simulation. We consider
several fixed values of $d$ ranging from $0.05$ to $0.6$ over roughly equally
spaced intervals. We fix the pilot sample size $m=5$ and the significance
level $\alpha=0.05.$ After implementing the sequential rule (9) with a
particular choice of $d,$ we obtain the confidence interval for $\lambda$ and
in turn report the interval for $\mu_{0}$, and finally an interval for the
average drying time of the droplet. Since the procedure has to be solved
analytically, all the analyses were carried out again using the
$``fitdistrplus"$ package in $R$.

Table 3. Analysis of the CA data using purely
---
sequential methodology (12) with $m=5,$ $\alpha=0.05$
$\underset{}{\overset{}{d}}$ | $N$ | CI for $\lambda$ | CI for $\mu_{0}$ | $\underset{}{\overset{}{\text{CI for Drying time (s)}}}$
$\overset{}{0.05}$ | $214$ | $(2.65,2.75)$ | $(0.34,0.36)$ | $(89.69,94.13)$
$0.1$ | $176$ | $(2.19,2.39)$ | $(0.39,0.42)$ | $(78.66,83.78)$
$0.2$ | $141$ | $(1.74,2.14)$ | $(0.43,0.51)$ | $(66.78,77.11)$
$0.3$ | $132$ | $(1.56,2.16)$ | $(0.43,0.56)$ | $(61.53,77.11)$
$0.4$ | $126$ | $(1.29,2.09)$ | $(0.44,0.65)$ | $(52.30,75.62)$
$0.5$ | $112$ | $(1.19,2.19)$ | $(0.42,0.69)$ | $(47.57,78.66)$
$0.6$ | $101$ | $(0.99,2.19)$ | $(0.42,0.78)$ | $(33.57,78.66)$

A few take away points from Table 3 are: as one goes on increasing $d$,
naturally, the width of the desired interval increases and as a result, less
number of observations are required to achieve a confidence level of $\alpha$
$(0.05$ in this case). Also, for increasing $N$, the intervals for the drying
time also increase and are seen to approach the actual estimated drying time
of $115.13$ seconds. But of course, a larger sample size comes with a cost and
hence one needs to strike a proper balance.

3\. Concluding remarks

In this paper we have established a sequential confidence interval methodology
to estimate the dynamic contact angle of a sessile saliva drop. This will help
the researchers and practitioners to build an idea about the growth of the
pandemic in general or by focusing on specific countries. Since a contact
angle has to be measured using a heavy-duty apparatus, a sequential rule also
appears to be beneficial by offering a reduction in time and cost. We
introduced a new circular model called as the half-circular wrapped-
exponential distribution to model the angles, which can only spread over half
a circle. This new model was seen to fit better than some of the existing ones
in the literature. Depending on the width $d$ of the interval fixed by the
experimenter, the mean contact angle of the droplet was seen to be between
0.41 and 0.56 radians or 23.49 and 32.08 degrees. On the other hand the drying
time of the saliva droplet was seen to be between 61 and 80 seconds.

References

[1] S. Balusamy, S. Banerjee and K. C. Sahu, Lifetime of sessile saliva
droplets in the context of SARS-CoV-2, Int. J. Heat Mass Trasf., 123, 105178
(2021).

[2] R. Mittal, R. Ni, and J.-H. Seo, The flow physics of COVID-19, J. Fluid
Mech. 894, F2 (2020).

[3] R. Bhardwaj and A. Agrawal, Likelihood of survival of coronavirus in a
respiratory droplet deposited on a solid surface, Phys. Fluids 32 (6), 061704
(2020).

[4] E. P. Vejerano and L.C. Marr, Physico-chemical characteristics of
evaporating respiratory fluid droplets, J. R. Soc. Interface 15 (139),
20170939 (2018).

[5] S. Chaudhuri, S. Basu, P. Kabi, V. R. Unni, and A. Saha, Modeling ambient
temperature and relative humidity sensitivity of respiratory droplets and
their role in determining growth rate of COVID-19 outbreaks, Phys. Fluids 32,
063309 (2020).

[6] E. Albert, B. Tegze, Z. Hajnal, D. Zámbó, D. P. Szekrényes, A. Déak, Z.
Hórvölgyi and N. Nagy, Robust contact angle determination for needle-in-drop
type measurements, ACS Omega 4, 18465-18471 (2019).

[7] K. V. Mardia and P. E. Jupp, Directional statistics, 2nd Ed. New York:
Wiley.

[8] J. S. Rao and A. Sengupta, Topics in circular statistics, New York: World
Scientific.

[9] A. V. D. Rao and S. V. S. Girija, Angular statistics, Boca Raton, CRC
Press.

[10] C. Ley and T. Verdebout, Applied directional statistics, Boca Raton, CRC
Press.

[11] A. Rambli, I. Mohamed, K. Shimizu and N. Ramli, A half-circular
distribution on a circle, Sains Malay., 48 (4), 887-892 (2019).

[12] J. S. Rao and T. J. Kozubowski, New families of wrapped distributions for
modeling skew circular data, Comm. in Stat. Theo. and Meth., 33 (9), 2059-2074
(2004).

[13] A. Yilmaz and C. Biçer, A new wrapped exponential distribution, Math.
Sci., 12, 285-293 (2018).

[14] S. Joshi and K. K. Jose, Wrapped Lindley distribution, Comm. in Stat.
Theo. and Meth., 47 (5) 1013-1021 (2018).

[15] S. Banerjee and N. Mukhopadhyay, A general sequential fixed-accuracy
confidence interval estimation methodology for a positive parameter:
Illustrations using health and safety data, Ann. of Inst. of Stat. Math., 68,
541-571 (2016).

[16] N. Mukhopadhyay and S. Banerjee, Purely sequential and two-stage bounded-
length confidence intervals for the Bernoulli parameter with illustrations
from health studies and ecology, P.K. Choudhary et al. (eds.), Ordered Data
Analysis, Modeling and Health Research Methods, Springer Proceedings in
Mathematics & Statistics 149.

[17] S. R. Bapat, On purely sequential estimation of an inverse Gaussian mean,
Metrika, 81, 1005-1024 (2018a).

[18] S. R. Bapat, Purely sequential fixed accuracy confidence intervals for
$P(X>Y)$ under bivariate exponential models, Am. J. of Math. and Manag. Sci.,
37, 386-400 (2018b).

[19] A. Khalifeh, E. Mahmoudi and A. Chaturvedi, Sequential fixed-accuracy
confidence intervals for the stress-strength reliability parameter for the
exponential distribution: two-stage sampling procedure, Comp. Stat.,
https://doi.org/10.1007/s00180-020-00957-5.

[20] S. Banerjee and N. Mukhopadhyay, A general sequential fixed-accuracy
confidence interval estimation methodology for a positive parameter:
illustrations using health and safety data, Ann. of the Inst. of Stat. Math.,
68, 541-570. (2016).