#  Distributed Iterative Hard Thresholding for Variable Selection in Tobit
Models

Changxin Yanga, Zhongyi Zhua and Heng Lianb
aDepartment of Statistics, Fudan University, Shanghai, China bDepartment of
Mathematics, City University of Hong Kong, Hong Kong, China

###### Abstract

While extensive research has been conducted on high-dimensional data and on
regression with left-censored responses, simultaneously addressing these
complexities remains challenging, with only a few proposed methods available.
In this paper, we utilize the Iterative Hard Thresholding (IHT) algorithm on
the Tobit model in such a setting. Theoretical analysis demonstrates that our
estimator converges with a near-optimal minimax rate. Additionally, we extend
the method to a distributed setting, requiring only a few rounds of
communication while retaining the estimation rate of the centralized version.
Simulation results show that the IHT algorithm for the Tobit model achieves
superior accuracy in predictions and subset selection, with the distributed
estimator closely matching that of the centralized estimator. When applied to
high-dimensional left-censored HIV viral load data, our method also exhibits
similar superiority.

Keywords— Censored regression; Distributed optimization; Hard thresholding;
High-dimension statistics; Linear convergence.

##  1 Introduction

The analysis of left-censored data is a significant statistical focus and has
attracted considerable research attention in recent years. It often arises due
to the lower detection limit of an assay, posing a shared challenge across
various fields such as biology, chemistry, and environmental sciences. For
example, biological assays used to measure Human Immunodeficiency Virus (HIV)
viral load in plasma may be limited in detecting concentrations below specific
thresholds. The presence of such missingness renders commonly used linear
regression methods ineffective. Furthermore, with advancements in modern data
collection, these challenges often manifest in high-dimensional scenarios. In
the context of HIV infection, there is a critical need to explore the
association between the number of viral loads and extremely high-dimensional
gene expression values. In addressing such challenges, various methodologies
have been proposed in prior research. Among these, the Tobit model has proved
to be valuable in modeling left-censored responses.

In recent years, significant progress has been made in the research on high
dimensional censored data. A useful strategy to tackle high-dimensionality
challenges involves constructing penalized estimators such as lasso-type
estimators. Müller and van de Geer (2016) and Zhou and Liu (2016) provide
theoretical insights into the least absolute deviation estimator with the
lasso penalty. Soret et al. (2018) propose the Lasso-regularized Buckley-James
least square algorithm, extending the estimator in Buckley and James (1979).
Jacobson and Zou (2024) were the first to consider the high-dimensional Tobit
model, which optimizes the likelihood function with a nonconvex penalty
(specifically with the SCAD penalty in Fan and Li (2001)). In this paper, we
concentrate on the IHT approach for variable selection. One advantage of using
IHT is that the user can directly specify the number of variables to be
retained, which may be useful in some scientific investigations.

IHT-style methods, which combine gradient descent with projection operations,
have gained popularity in the literature for sparse recovery. Various
algorithms have been proposed, such as standard IHT introduced in Blumensath
and Davies (2009), GraDeS in Garg and Khandekar (2009), and Hard Thresholding
Pursuit (HTP) in Foucart (2011). Jain et al. (2014) demonstrated that the IHT
procedure can achieve linear convergence to the optimal solution under the
conditions of strong convexity and strong smoothness in high-dimensional
settings. Extending this result, Wang et al. (2023) considered nonsmooth loss
functions under the less restrictive assumption of a locally positive-definite
population Hessian. We will demonstrate that the linear convergence result
also holds for the high-dimensional Tobit model.

Additionally, we consider a scenario where data are distributed across
multiple locations and develop a distributed variant of our method. With
advancements in data collection technology, gathering similar types of data
from various regions has become increasingly common. Typically, directly
aggregating all data at a central site faces practical hurdles related to
storage, communication, and privacy concerns. To tackle these challenges,
divide-and-conquer approaches are commonly employed. Many early approaches
utilize one-shot methods, where estimators computed on local machines are
transmitted to a central node and aggregated to form a global estimate, as in
Zhang et al. (2013), Lee et al. (2017). However, such methods suffer from some
drawbacks, as discussed by Jordan et al. (2019), who proposed a communication-
efficient distributed algorithm. Fortunately, the gradient descent-type update
in the IHT method is conducive to integration with this algorithm. Building
upon this concept, we introduce a communication-efficient estimator with an
IHT-type update process. To our knowledge, this represents the first
exploration of the high-dimensional Tobit model in a distributed setting. Our
theoretical analysis and numerical results demonstrate that this estimator’s
convergence rate aligns with that achieved when pooling all data together
while incurring low communication costs.

This article is structured as follows. In Section 2, we review the Tobit model
and then develop both local (centralized) IHT and the distributed IHT. Section
3 is dedicated to our theoretical analysis. Section 4 draws some conclusions.
All proofs are relegated to the Appendix.

##  2 The Tobit model and IHT

###  2.1 Local IHT for Tobit model

We consider the Tobit model, introduced in Tobin (1958), a crucial tool for
modeling a left-censored response. We assume the existence of a latent
response variable $y^{*}$ such that $y=\max\left\\{y^{*},c_{0}\right\\}$,
where $y^{*}$ follows a linear model

$y^{*}=\mathbf{x}^{\prime}\bm{\beta}^{*}+\epsilon.$

Here,
$\mathbf{x}=\left(1,x_{1},\ldots,x_{d}\right)^{\prime}\in\mathbb{R}^{d+1}$,
$\bm{\beta}^{*}=\left(\beta_{0}^{*},\beta_{1}^{*},\ldots,\beta_{d}^{*}\right)^{\prime}\in\mathbb{R}^{d+1}$,
and $\epsilon\sim\mathcal{N}\left(0,(\sigma^{*})^{2}\right)$. We focus on a
high-dimensional scenario where the dimension $d$ is much larger than the
sample size $n$, with $\left\|\bm{\beta}^{*}\right\|_{0}=s_{0}$. Without loss
of generality, we assume $c_{0}=0$ throughout the following discussions.

Let $\\{\mathbf{z}_{i}=\left(\mathbf{x}_{i},y_{i}\right)\\}_{i=1}^{n}$
represent an independent and identically distributed (i.i.d.) sample from the
Tobit model, and define the indicator $d_{i}={I}_{y_{i}>0}$. The likelihood
function for the censored response in the Tobit model can be expressed as

$\displaystyle
L_{n}\left(\bm{\beta},\sigma^{2},\\{\mathbf{z}_{i}\\}_{i=1}^{n}\right)=\prod_{i=1}^{n}\left[{1\over\sqrt{2\pi}\sigma}\exp\left\\{-{1\over
2\sigma^{2}}\left(y_{i}-\mathbf{x}_{i}^{\prime}\bm{\beta}\right)^{2}\right\\}\right]^{d_{i}}\left[\Phi\left({-\mathbf{x}_{i}^{\prime}\bm{\beta}\over\sigma}\right)\right]^{1-d_{i}},$

where $\Phi(\cdot)$ is the standard normal cumulative distribution function.
While ${L}_{n}\left(\bm{\beta},\sigma^{2}\right)$ is not concave in
$\left(\bm{\beta},\sigma^{2}\right)$, Olsen (1978) found that the
reparameterization $\bm{\delta}=\mbox{\boldmath$\beta$}/\sigma$ and
$\gamma^{2}=\sigma^{-2}$ results in a concave log-likelihood. Denoting
$\bm{\theta}=(\bm{\delta},\gamma)$, after dropping ignorable constants, the
negative log-likelihood is given by

$\displaystyle\widehat{\mathcal{L}}(\bm{\theta})$ $\displaystyle={1\over
n}\sum_{i=1}^{n}{\mathcal{L}}(\bm{\theta},\mathbf{z}_{i})$
$\displaystyle={1\over n}\sum_{i=1}^{n}d_{i}\left[-{\rm log}(\gamma)+{1\over
2}\left(\gamma
y_{i}-\mathbf{x}_{i}^{\prime}\bm{\delta}\right)^{2}\right]-\left(1-d_{i}\right){\rm
log}\left(\Phi\left(-\mathbf{x}_{i}^{\prime}\bm{\delta}\right)\right).$

In the high-dimensional setting that $d$ is large, and maybe even larger than
$n$, variable selection is desired. In this paper, we consider the hard
constraint $\|\mbox{\boldmath$\delta$}\|_{0}\leq s$ where $s$ is a known upper
bound of $s_{0}$. To optimize this constrained objective function, we use the
IHT method, which is designed for high-dimensional regression. This method
combines gradient descent with projection operations, making it
computationally efficient. Furthermore, we slightly modify it for the Tobit
model to prevent $\gamma$ from getting too small or even negative.
Specifically, denoting the initial value as $\bm{\theta}^{0}$, for
$t=0,1,\ldots$, our IHT algorithm can be formulated as

$\bm{\theta}^{t+1}=P_{s,C^{*}}\left(\bm{\theta}^{t}-\eta\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{t}\right)\right)\doteq\left\\{\begin{array}[]{l}\bm{\delta}^{t+1}=P_{s}\left(\bm{\delta}^{t}-\eta\nabla_{\bm{\delta}}\widehat{\mathcal{L}}\left(\bm{\theta}^{t}\right)\right)\\\
\gamma^{t+1}=\mathcal{T}_{C^{*}}(\gamma^{t}-\eta\nabla_{\gamma}\widehat{\mathcal{L}}\left(\bm{\theta}^{t}\right))\end{array}\right.,$
(1)

where $\eta>0$ is the step size, $P_{s}(\mathbf{\cdot})$ is the projection
operator that retains only $s$ entries with the largest absolute values while
setting other entries to zero, and

$\mathcal{T}_{C^{*}}(\gamma^{t})=\begin{cases}\gamma^{t}-\eta\nabla_{\gamma}\widehat{\mathcal{L}}\left(\bm{\theta}^{t}\right)&\text{if
}\gamma^{t}-\eta\nabla_{\gamma}\widehat{\mathcal{L}}\left(\bm{\theta}^{t}\right)\geq
C^{*}\\\ C^{*}&\text{otherwise}\end{cases}.$

Algorithm 1 Local IHT for Tobit Regression

1:sparsity level $s$, number of iterations $T$, step size $\eta$, lower bound
$C^{*}$.

2:Initialize $\bm{\theta}^{0}_{local}$

3:for $t=0,\ldots,T-1$ do

4: Calculate the gradient
$\nabla\widehat{\mathcal{L}}(\bm{\theta}^{t}_{local})$

5: Update $\bm{\theta}^{t+1}_{local}$ as
$\bm{\theta}^{t+1}_{local}=P_{s,C^{*}}\left(\bm{\theta}^{t}_{local}-\eta\nabla\widehat{\mathcal{L}}\left(\bm{\theta}_{local}^{t}\right)\right)$

6:end for

7:return
$\bm{\theta}^{T}_{local}=(\bm{\delta}^{T}_{local},\gamma^{T}_{local})$

Note that we have excluded $\gamma$ from the projection operation, as $\gamma$
is not part of the parameter for variable selection. Additionally, we
introduce a positive lower bound denoted as $C^{*}$ to prevent $\gamma\leq 0$
during the update process.

In practical applications, the values of $s$ and $C^{*}$ are typically unknown
and are treated as tuning parameters. There are two main approaches for
determining $s$. First, the user can directly specify it based on a predefined
desired sparsity, such as selecting the top 10 variables of interest. An
alternative is to treat sparsity as a tuning parameter that can be optimized
through cross-validation. $C^{*}$ can be set to a sufficiently small value to
ensure it is smaller than $\gamma^{*}$. We present the entire IHT algorithm in
Algorithm 1 (called local IHT to distinguish it from the distributed cases).
Theorem 1 in Section 3 provides a bound on the estimator error for the local
IHT algorithm (after a sufficient number of iterations) as
$O\left(\sqrt{{s{\rm log}(n\vee d)\over n}}\right)$, which is well-known as
the minimax rate for linear regression under the sparsity condition (if
$s_{0}\ll d$, and $s$ and $s_{0}$ are of the same order).

###  2.2  Distributed IHT for Tobit model

In this section, we delve into the distributed implementation of the high-
dimensional Tobit model across $M$ machines. This scenario holds significance
in settings characterized by large volumes of training data. In such contexts,
aggregating all raw data directly becomes impractical due to constraints like
limited storage capacity, high communication costs, and privacy
considerations. Hence, we are motivated to explore the development of a
communication-efficient distributed IHT method.

For simplicity, we assume that each machine stores data of the same sample
size $n$. Let $\mathcal{H}_{m}$ denote the subsample with a sample size of $n$
stored in the $m$-th machine $\mathcal{M}_{m}$, for $m=1,\ldots,M$, and denote
the global sample size $N=nM$. The local and global negative log-likelihood
functions take the following form:

$\begin{gathered}\mathcal{L}_{m}(\bm{\theta})={1\over
n}\sum_{i\in\mathcal{H}_{m}}{\mathcal{L}}(\bm{\theta},\mathbf{z}_{i}),\\\
\mathcal{L}_{N}(\bm{\theta})={1\over
M}\sum_{m=1}^{M}\mathcal{L}_{m}(\bm{\theta}).\end{gathered}$

Our target is to find

${\bm{\theta}}^{*}=\underset{\bm{\theta}:\|\bm{\theta}\|_{0}\leq
s+1}{\operatorname{argmin}}\mathcal{L}^{*}(\bm{\theta})\doteq
E(\mathcal{L}_{N}(\bm{\theta})).$ (2)

When direct aggregation of all data to optimize the global loss function is
infeasible, a common approach is to approximate the global loss function via
inter-machine communication. Subsequently, this approximation is used for
optimization based on local data. Jordan et al. (2019) demonstrated that this
approach could achieve optimal statistical precision with just a few rounds of
communication. Following suit, we replace the original global negative log-
likelihood function with an approximation. Specifically, given the current
estimate $\overline{\bm{\theta}}$, we first collect the gradients
$\nabla\mathcal{L}_{m}(\overline{\bm{\theta}})$ from different machines and
then approximate the global negative log-likelihood function as

$\widetilde{\mathcal{L}}(\bm{\theta}):=\mathcal{L}_{1}(\bm{\theta})-\left\langle\bm{\theta},\nabla\mathcal{L}_{1}(\overline{\bm{\theta}})-\nabla\mathcal{L}_{N}(\overline{\bm{\theta}})\right\rangle,$
(3)

where $\mathcal{L}_{1}(\bm{\theta})$ is the local empirical negative log-
likelihood function on the first machine $\mathcal{M}_{1}$ (treated as the
central machine that can communicate with all others) and
$\nabla\mathcal{L}_{N}(\overline{\bm{\theta}})={1\over
M}\sum_{i=1}^{M}\nabla\mathcal{L}_{m}(\overline{\bm{\theta}})$. The global
gradient for the approximate loss is
$\nabla\widetilde{\mathcal{L}}(\bm{\theta}):=\nabla\mathcal{L}_{1}(\bm{\theta})-(\nabla\mathcal{L}_{1}(\overline{\bm{\theta}})-\nabla\mathcal{L}_{N}(\overline{\bm{\theta}})).$
Thus, based on the approximate global negative log-likelihood function, the
IHT algorithm is, starting from
$\mbox{\boldmath$\theta$}^{0}=\overline{}\mbox{\boldmath$\theta$}$,

$\bm{\theta}^{t+1}=P_{s,C^{*}}\left(\bm{\theta}^{t}-\eta\nabla\widetilde{\mathcal{L}}(\bm{\theta}^{t})\right).$
(4)

After (4) converges, ${\cal M}_{1}$ sends the new estimate to all other
machines to obtain the updated approximate global negative log-likelihood
function. Algorithm 2 presents the details of the distributed IHT algorithm
for the Tobit model. Our theoretical results will show that the distributed
IHT method can achieve the same statistical rate as the centralized version.
Our simulation results further validate that it offers comparable statistical
precision compared to its non-distributed counterpart.

Algorithm 2 Distributed IHT for Tobit Regression

1:sparsity level $s$, number of iterations in the inner loop $T$, number of
iterations for the outer loop $Q$, step size $\eta$, lower bound $C^{*}$.

2:Initialize $\bm{\theta}^{0}_{dis}$.

3:for $q=0,\ldots,Q-1$ do

4: Let $\bar{\bm{\theta}}=\bm{\theta}^{q}_{dis}$ and ${\cal M}_{1}$ sends
$\bar{\bm{\theta}}$ to all other machines.

5: Each machine $\mathcal{M}_{m}$ computes
$\nabla\mathcal{L}_{m}(\bar{\bm{\theta}})$ and send it to ${\cal M}_{1}$.

6: $\mathcal{M}_{1}$ calculates
$\nabla\mathcal{L}_{N}(\bar{\bm{\theta}})={1\over
M}\sum_{m=1}^{M}\nabla\mathcal{L}_{m}(\bar{\bm{\theta}}).$

7: On ${\cal M}_{1}$, let ${\bm{\theta}}^{0}_{inner}=\bar{\bm{\theta}}$ and
execute the following inner loop.

8: for $t=0,\ldots,T-1$ do
$\bm{\theta}^{t+1}_{inner}=P_{s,C^{*}}\left(\bm{\theta}_{inner}^{t}-\eta\nabla\widetilde{\mathcal{L}}(\bm{\theta}_{inner}^{t})\right).$

9: end for

10: Let $\bm{\theta}^{q+1}_{dis}=\bm{\theta}^{T}_{inner}$

11:end for

12:return $\bm{\theta}^{Q}_{dis}=(\bm{\delta}^{Q}_{dis},\gamma^{Q}_{dis})$

##  3 Convergence and statistical rates

We present convergence guarantees for both the local IHT and the distributed
IHT in this section. In the rest of the paper, $c$ and $C$ denote positive
constants independent of $n$, $d$, $M$, whose values may vary across
instances. In the following, for positive sequences $a_{n}$ and $b_{n}$,
$a_{n}=O\left(b_{n}\right)$ or $a_{n}\lesssim b_{n}$ means $a_{n}/b_{n}\leq c$
for some constant $c$, $a_{n}=o\left(b_{n}\right)$ means
$a_{n}/b_{n}\rightarrow 0$, $a_{n}\asymp b_{n}$ means $a_{n}=O(b_{n})$ and
$b_{n}=O(a_{n})$, $a_{n}\ll b_{n}$ or $b_{n}\gg a_{n}$ is the same as
$a_{n}=o\left(b_{n}\right)$. We make the following assumptions.

  * $\bullet$

Assumption 1 (Truncation and signal strength): $C^{*}$ is a constant that
satisfies $\left.0<C^{*}<\gamma^{*}\right.$.
$\|\mbox{\boldmath$\beta$}^{*}\|_{0}=s_{0}$ and
$\|\mbox{\boldmath$\beta$}^{*}\|_{2}\leq C$.

  * $\bullet$

Assumption 2 (Sparse eigenvalue condition): There exist positive constants
$C_{1}$ and $\tilde{C}_{2}$ such that:

    1. (i)

$\displaystyle C_{1}\leq
E\left[\mathbf{a}^{\prime}\begin{pmatrix}I_{y>0}\mathbf{x}\mathbf{x}^{\prime}&-I_{y>0}\mathbf{x}y\\\
-I_{y>0}y\mathbf{x}^{\prime}&I_{y>0}y^{2}\end{pmatrix}\mathbf{a}\right]<E\left[\mathbf{a}^{\prime}\begin{pmatrix}\mathbf{x}\mathbf{x}^{\prime}&-I_{y>0}\mathbf{x}y\\\
-I_{y>0}y\mathbf{x}^{\prime}&I_{y>0}y^{2}\end{pmatrix}\mathbf{a}\right]\leq\tilde{C}_{2},$

for all $(2s+s_{0}+1)$-sparse unit vectors $\mathbf{a}$.

    2. (ii)

$\left(1+\sqrt{{s_{0}\over s}}\right){C_{2}-C_{1}\over C_{2}+C_{1}}<1$ holds,
where $C_{2}\dot{=}\tilde{C}_{2}+{1\over\left(C^{*}\right)^{2}}$, and
$s=O(s_{0})$.

  * $\bullet$

Assumption 3 (Sub-Gaussian designs): For some constant $C$, $\mathbf{x}_{i}$
are iid and $C$-sub-Gaussian, denoted as $\operatorname{subG}(C)$. That is,
for all $\alpha\in\mathbb{R}^{d}$,

$\mathbb{E}[\exp(\alpha^{\prime}\mathbf{x})]\leq\exp\left({C^{2}\|\alpha\|_{2}^{2}\over
2}\right).$

  * $\bullet$

Assumption 4 (Proper initialization): $\|\bm{\theta}_{int}\|_{0}\leq(s+1)$ and
$\gamma_{int}\geq C^{*}$.

Remark 1: To satisfy Assumption 1, we need to use a sufficiently small
$C^{*}$. The sparse eigenvalue assumption is frequently used in high-
dimensional data analysis. Note that here the assumption is imposed on the
population quantity so even the standard eigenvalue assumption without
constraining ${\bf a}$ to be sparse is reasonable. Examining the form of the
population Hessian matrix $H(\bm{\theta})$ (seen in the Appendix), it is easy
to see that Assumptions 1 and 2 (i) actually imply the Hessian matrix has
eigenvalues between $C_{1}$ and $C_{2}$. Assumption 3 is also a standard
assumption in high-dimensional regression analysis, making it possible to
derive tight bounds on different quantities (Jordan et al., 2019; Belloni and
Chernozhukov, 2011). Combined with the assumption
$\|\mbox{\boldmath$\beta$}^{*}\|_{2}\leq C$ and the fact that the noise is
Gaussian, it also implies $y_{i}$ is sub-Gaussian. Finally, a trivial initial
value satisfying assumption 4 is
$\mbox{\boldmath$\theta$}_{int}=(\mathbf{0},C^{*})$.

The following Theorem provides a bound for the $\ell_{2}$ estimation error of
the local IHT algorithm, confirming the linear convergence of the IHT
algorithm for the Tobit model.

######  Theorem 1.

(Local IHT for Tobit) Under Assumptions 1–4, choosing the step size
$\eta={2\over C_{1}+C_{2}}$,

$\displaystyle\|\mbox{\boldmath$\theta$}_{local}^{t+1}-\mbox{\boldmath$\theta$}^{*}\|_{2}$
$\displaystyle\leq$ $\displaystyle\left(\left(1+\sqrt{{s_{0}\over
s}}\right){C_{2}-C_{1}\over C_{2}+C_{1}}+Cs\sqrt{{{\rm log}(d\vee n)\over
n}}\right)\|\mbox{\boldmath$\theta$}_{local}^{t}-\mbox{\boldmath$\theta$}^{*}\|_{2}$
$\displaystyle+C\sqrt{s}\|\mbox{\boldmath$\theta$}_{local}^{t}-\mbox{\boldmath$\theta$}^{*}\|_{2}^{2}+C\sqrt{{s{\rm
log}(d\vee n)\over n}},$

with probability at least $1-(d\vee n)^{-C}$. In particular, if ${s^{2}{\rm
log}(d\vee n)\over n}=o(1)$ and
$\sqrt{s}\|\mbox{\boldmath$\theta$}_{\text{int}}-\mbox{\boldmath$\theta$}^{*}\|=o(1)$,
after $T_{1}\asymp{\rm
log}\left({n\left\|\bm{\theta}_{\text{int}}-\bm{\theta}^{*}\right\|^{2}_{2}\over
s{\rm log}(d\vee n)}\right)$ iterations, we have

$\left\|\bm{\theta}_{\text{local}}^{T_{1}}-\bm{\theta}^{*}\right\|_{2}\lesssim\sqrt{{s{\rm
log}(d\vee n)\over n}},$

with probability at least $1-(d\vee n)^{-C}$.

The specific value chosen for $\eta$ is not the only option, and we opt for it
for convenience in the proof. Our proof reveals that the requirement for
$\eta$ and $s$ is actually that $\left(1+\sqrt{{s_{0}\over
s}}\right)\max\left\\{\left|1-\eta C_{1}\right|,\left|1-\eta
C_{2}\right|\right\\}<1$. For initialization, we recommend the Tobit-lasso
estimator as it offers computational simplicity and has been shown to be
consistent in Jacobson and Zou (2024).

Theorem 1 establishes the convergece of the local IHT algorithm. The following
Theorem concerns the convergence of the inner loop in the distributed setting,
which is very similar to Theorem 1. This is not surprising since the
approximate likelihood is expected to be close to the true likelihood.

######  Theorem 2.

(Inner loop for Distributed IHT) Given $\bar{\bm{\theta}}=\bm{\theta}_{int}$
in Algorithm 2, under Assumptions 1–4, choosing the step size $\eta={2\over
C_{1}+C_{2}}$, if ${s^{2}{\rm log}(d\vee n)\over n}=o(1)$ and
$\sqrt{s}\|\bar{}\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{*}\|=o(1)$,
then after $T_{2}\asymp{\rm
log}\left({n\|\bar{}\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{*}\|^{2}\over
s{\rm log}(d\vee n)}\right)$ iterations, we have

$\left\|\bm{\theta}^{T_{2}}_{\text{inner}}-\bm{\theta}^{*}\right\|_{2}\lesssim\sqrt{s}\left(\sqrt{{s{\rm
log}(d\vee n)\over
n}}\left\|\bar{\bm{\theta}}-\bm{\theta}^{*}\right\|_{2}+\|\bar{\bm{\theta}}-\bm{\theta}^{*}\|_{2}^{2}+\sqrt{{{\rm
log}(d\vee n)\over N}}\right),$ (5)

with probability at least $1-(d\vee n)^{-C}$.

The above result shows that after the inner loop ends, we have possibly a
better estimate compared to when the loop starts. Combining such bounds over
multiple stages $Q$, we get the following result.

######  Theorem 3.

(Distributed IHT for Tobit) Under the assumptions of Theorem 2, after
$Q\asymp{{\rm log}N\over{\rm log}n}$ iterations in the outer loop (with the
number of iterations in the inner loop as in the last statement of Theorem 2),
with probability at least $1-(d\vee n)^{-C}$, we have

$\left\|\bm{\theta}^{Q}_{dis}-\bm{\theta}^{*}\right\|_{2}\lesssim\sqrt{{s{\rm
log}(d\vee n)\over N}}.$

##  4 Conclusion

In this paper, we delved into applying the IHT algorithm to the high-
dimensional Tobit model and expanded its utility to distributed environments
spanning multiple machines. This marks the first effort to study the Tobit
model under the distributed setting in high-dimensional scenarios. We have
shown that achieving a minimax convergence rate for local and distributed
estimators is feasible. Simulation results validate our theoretical findings
and demonstrate the good performance of our proposed method across diverse
settings. When applied to high-dimensional left-censored HIV viral load data,
the IHT method also effectively produced accurate predictions.

Regarding future directions, it is possible to explore extending the IHT
algorithm to decentralized distributed algorithms, such as distributed
gradient descent (DGD) as introduced in Nedic and Ozdaglar (2009). DGD
algorithms are extensively studied in the literature and offer advantages such
as avoiding a single point of failure, addressing the limitations of
centralized frameworks. Relaxing the requirement for sufficient local sample
size by leveraging information from neighboring nodes in decentralized
networks is also possible, as discussed in Ji et al. (2023). Exploring the
combination with other distributed methods also remains an intriguing area for
future research. Moreover, we can also consider a more complex class of
missing data models known as Zero-Inflated models, where missing values are
characterized as arising from a two-component mixture model. The presence of
latent components in missing values requires a more careful analysis.

## Appendix: Proofs

The proof strategies for Theorem 1 and Theorem 2 share similarities.
Therefore, we will first provide a detailed proof for Theorem 2, which is
expected to be more challenging. We introduce some additional notations used
in the proof.

  * $\bullet$

Let $\phi(\cdot)$ and $\Phi(\cdot)$ denote the standard normal probability
density function (PDF) and cumulative distribution function (CDF),
respectively. Let $g(a)={\phi(a)\over\Phi(a)}$. Define
$h(a)=g(a)(a+g(a))=-g^{\prime}(a)$. According to Sampford (1953) and Jacobson
and Zou (2024), we have $0<h(a)<1$ and $-4.3<h^{\prime}(a)<0$.

  * $\bullet$

Let
$\nabla\mathcal{L}_{1}^{t}=\nabla\mathcal{L}_{1}\left(\bm{\theta}^{t}\right)$
denote the negative log-likelihood of the first machine, and similarly for
$\nabla\tilde{\cal L}^{t}=\nabla\tilde{\mathcal{L}}(\bm{\theta}^{t})$. The
gradient of the log-likelihood $\nabla\mathcal{L}_{1}(\bm{\theta})$ can be
expressed as

${1\over n}\sum_{i\in{\cal
H}_{1}}\left[\begin{array}[]{c}-I_{y_{i}>0}\mathbf{x}_{i}\left(\gamma
y_{i}-\mathbf{x}_{i}^{\prime}\mbox{\boldmath$\delta$}\right)+I_{y_{i}\leq
0}\left(\mathbf{x}_{i}^{\prime}g\left(-\mathbf{x}^{\prime}_{i}\bm{\delta}\right)\right)\\\
-I_{y_{i}>0}\left(\gamma^{-1}-y_{i}^{\prime}\left(\gamma
y_{i}-\mathbf{x}_{i}^{\prime}\mbox{\boldmath$\delta$}\right)\right)\end{array}\right].$
(6)

  * $\bullet$

Let $H(\bm{\theta})$ denote the population Hessian matrix, given by

$E\left(\begin{bmatrix}I_{y>0}\mathbf{xx}^{\prime}+I_{y\leq
0}\mathbf{xx}^{\prime}h\left(-\mathbf{x}^{\prime}\bm{\delta}\right)&-I_{y>0}\mathbf{x}y\\\
-I_{y>0}y\mathbf{x}^{\prime}&I_{y>0}(y^{2}+\gamma^{-2})\
\end{bmatrix}\right).$ (7)

  * $\bullet$

Define $I^{t}=S^{t}\cup S_{0}$, where $S^{t}$ and $S_{0}$ represent the
support (indices of nonzero entries) of $\bm{\theta}^{t}$ and the true
$\bm{\theta}^{*}$, respectively.

  * $\bullet$

Given a vector $\bm{\theta}$, a matrix $B$, and an index set
$I\subseteq\\{1,\ldots,d+2\\}$, $\bm{\theta}_{I}$ denotes the vector obtained
from $\bm{\theta}$ by setting components not in $I$ to zero, while $B_{I}$
represents the matrix obtained from $B$ by setting the components of rows not
in $I$ to zero. $|I|$ represents the cardinality of the set $I$.

  * $\bullet$

We say that a random variable $X$ is sub-exponential with parameters $C$ if
the $L^{p}$ norm of $X$ satisfy

$\|X\|_{L^{p}}=\left({E}|X|^{p}\right)^{1/p}\leq Cp,\quad\text{ for all }p\geq
1.$

  * $\bullet$

Define
$\mathbf{x}_{\max}=\max_{i\in\\{1,\dots,n\\}}||\mathbf{x}_{i}||_{\infty}$, and
$y_{\max}=\max_{i\in\\{1,\dots,n\\}}|y_{i}|$.

### Proof of Theorem 2

In the inner loop, we have

$\displaystyle e_{t+1}$
$\displaystyle\doteq\left\|\bm{\theta}^{t+1}-\bm{\theta}^{*}\right\|_{2}$
$\displaystyle=\left\|\left(P_{s,C^{*}}\left(\bm{\theta}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}$
$\displaystyle\leq\left\|\left(P_{s}\left(\bm{\theta}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)\right)_{I^{t+1}}-\left(\bm{\theta}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\left\|\left(\bm{\theta}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}$
$\displaystyle\leq\left(1+\sqrt{{\left|I^{t+1}\right|-s\over\left|I^{t+1}\right|-s_{0}}}\right)\left\|\left(\bm{\theta}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}.$

The first inequality above holds because, under Assumption 1, if
$\gamma^{t}-\eta\nabla_{\gamma}\tilde{\mathcal{L}}\left(\bm{\theta}^{t}\right)<C^{*}$,
then it is evident that
$\left|\gamma^{*}-C^{*}\right|\leq\left|\gamma^{*}-\gamma^{t}-\eta\nabla_{\gamma}\tilde{\mathcal{L}}\left(\bm{\theta}^{t}\right)\right|$,
meaning removing the thresholding with $C^{*}$ can only increase the error.
The 2nd inequality is due to Lemma 1 in Jain et al. (2014).

Then, by adding and subtracting terms, we have

$\displaystyle\left\|\left({\bm{\theta}}^{t}-\eta\nabla\tilde{\mathcal{L}}^{t}\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}$
$\displaystyle\leq\left\|{\bm{\theta}}_{I^{t+1}}^{t}-\bm{\theta}^{*}_{I^{t+1}}-\eta
E\left(\nabla\tilde{\mathcal{L}}^{t}-\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|E\left(\nabla\tilde{\mathcal{L}}^{t}-\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}-\left(\nabla\tilde{\mathcal{L}}^{t}-\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|\left(\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle=\left\|{\bm{\theta}}_{I^{t+1}}^{t}-\bm{\theta}^{*}_{I^{t+1}}-\eta
E\left(\nabla\mathcal{L}_{1}^{t}-\nabla{\mathcal{L}_{1}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|E\left(\nabla\mathcal{L}_{1}^{t}-\nabla{\mathcal{L}_{1}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}-\left(\nabla\mathcal{L}_{1}^{t}-\nabla{\mathcal{L}_{1}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|\left(\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\doteq||P_{1}||_{2}+\eta||P_{2}||_{2}+\eta||P_{3}||_{2}.$

For $P_{1}$, by Assumption 2 and Taylor’s expansion, we get

$\displaystyle\left\|\bm{\theta}_{I^{t+1}}^{t}-\bm{\theta}^{*}_{I^{t+1}}-\eta\left(E\nabla\mathcal{L}_{1}^{t}-E\nabla\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
(8) $\displaystyle\leq\|(I-\eta H(\bm{\alpha}))_{I^{t+1}\cup
S^{t}}\left(\bm{\theta}_{I^{t+1}\cup S^{t}}^{t}-\bm{\theta}^{*}_{I^{t+1}\cup
S^{t}}\right)\|_{2}$ $\displaystyle\leq\max\left\\{\left|1-\eta
C_{1}\right|,\left|1-\eta
C_{2}\right|\right\\}\left\|\bm{\theta}^{t}-\bm{\theta}^{*}\right\|_{2},$

where $\bm{\alpha}$ lies on the line between $\bm{\theta}^{t}$ and
$\bm{\theta}^{*}$.

For $P_{2}$, by Lemma 2, with probability at least $1-(d\vee n)^{-C}$, we have

$\displaystyle\|P_{2}\|_{\infty}$ $\displaystyle\lesssim\sqrt{{s{\rm
log}(d\vee n)\over n}}e_{t}+e_{t}^{2},$ (9)

and thus $\|P_{2}\|_{2}\lesssim s\sqrt{{{\rm log}(d\vee n)\over
n}}e_{t}+\sqrt{s}e_{t}^{2}$.

For $P_{3}$, by Lemma 3, we have, with probability at least $1-(d\vee
n)^{-C}$,

$\displaystyle\|P_{3}\|_{2}$ $\displaystyle\lesssim s\sqrt{{{\rm log}(d\vee
n)\over n}}e_{0}+\sqrt{s}e_{0}^{2}+\sqrt{{s{\rm log}(d\vee n)\over N}}.$

Therefore, by choosing $\eta={2\over C_{2}+C_{1}}$ and combining bounds for
$P_{1}$, $P_{2}$, and $P_{3}$, we have

$\displaystyle e_{t+1}$ $\displaystyle\leq$
$\displaystyle\left(\left(1+\sqrt{{s_{0}\over s}}\right){C_{2}-C_{1}\over
C_{2}+C_{1}}+Cs\sqrt{{{\rm log}(d\vee n)\over
n}}\right)e_{t}+\sqrt{s}e_{t}^{2}$ $\displaystyle+Cs\sqrt{{{\rm log}(d\vee
n)\over n}}e_{0}+C\sqrt{s}e_{0}^{2}+\sqrt{{s{\rm log}(d\vee n)\over N}},$

with probability at least $1-(d\vee n)^{-C}$.

Thus, if $s\sqrt{{{\rm log}(d\vee n)\over n}}=o(1)$ and $\sqrt{s}e_{0}=o(1)$,
after $T\asymp\left({\rm log}\left({ne_{0}\over s{\rm log}(d\vee
n)}\right)\right)$ iterations, with probability at least $1-(d\vee n)^{-C}$,

$e_{T}\lesssim\sqrt{s}\left(\sqrt{{s{\rm log}(d\vee n)\over
n}}\left\|\bar{\bm{\theta}}-\bm{\theta}^{*}\right\|_{2}+\|\bar{\bm{\theta}}-\bm{\theta}\|_{2}^{2}+\sqrt{{{\rm
log}(d\vee n)\over N}}\right).$ (10)

The proof is completed.

### Proof of Theorem 1

In fact, the proof of Theorem 1 is essentially the same as the proof of
Theorem 2. Our goal is to provide an upper bound for

$\displaystyle e_{t+1}$
$\displaystyle\doteq\left\|\bm{\theta}^{t+1}-\bm{\theta}^{*}\right\|_{2}\leq\left(1+\sqrt{{\left|I^{t+1}\right|-s\over\left|I^{t+1}\right|-s_{0}}}\right)\left\|\left(\bm{\theta}^{t}-\eta\nabla\widehat{\mathcal{L}}^{t}\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}.$

We have

$\displaystyle\left\|\left({\bm{\theta}}^{t}-\eta\nabla\widehat{\mathcal{L}}^{t}\right)_{I^{t+1}}-\bm{\theta}^{*}_{I^{t+1}}\right\|_{2}$
$\displaystyle\leq\left\|{\bm{\theta}}_{I^{t+1}}^{t}-\bm{\theta}^{*}_{I^{t+1}}-\eta
E\left(\nabla\widehat{\mathcal{L}}^{t}-\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|E\left(\nabla\widehat{\mathcal{L}}^{t}-\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}-\left(\nabla\widehat{\mathcal{L}}^{t}-\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\quad+\eta\left\|\left(\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right)_{I^{t+1}}\right\|_{2}$
$\displaystyle\doteq||P_{1}||_{2}+\eta||P_{2}||_{2}+\eta||P_{3}||_{2}.$

For $P_{1}$, by Assumption 2 and Taylor’s expansion, we can obtain that

$\left\|P_{1}\right\|_{2}\leq\max\left\\{\left|1-\eta
C_{1}\right|,\left|1-\eta C_{2}\right|\right\\}e_{t}.$

For $P_{2}$, similar to (9), with probability at least $1-(d\vee n)^{-C}$, we
have

$\displaystyle\left\|P_{2}\right\|_{2}$ $\displaystyle\lesssim s\sqrt{{{\rm
log}(d\vee n)\over n}}e_{t}+\sqrt{s}e_{t}^{2}.$

For $P_{3}$, as we have shown in the proof of Lemma 3,
$\left\|\nabla\widehat{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right\|_{2}\lesssim\sqrt{{s{\rm
log}(d\vee n)\over n}}$ with probability at least $1-(d\vee n)^{-C}$.

Combining the bounds and choosing $\eta={2\over C_{2}+C_{1}}$, we get

$\displaystyle e_{t+1}$ $\displaystyle\leq\left(\left(1+\sqrt{{s_{0}\over
s}}\right){C_{2}-C_{1}\over C_{2}+C_{1}}+Cs\sqrt{{{\rm log}(d\vee n)\over
n}}\right)e_{t}+C\sqrt{s}e_{t}^{2}+C\sqrt{{s{\rm log}(d\vee n)\over n}},$

with probability of at least $1-(d\vee n)^{-C}$.

Thus if $s\sqrt{{{\rm log}(d\vee n)\over n}}=o(1)$ and $\sqrt{s}e_{0}=o(1)$,
after $T_{1}\asymp{\rm
log}\left({n\left\|\bm{\theta}^{0}-\bm{\theta}^{*}\right\|^{2}_{2}\over s{\rm
log}(d\vee n)}\right)$ iterations , with probability at least $1-(d\vee
n)^{-C}$, we have

$e_{T}\lesssim\sqrt{{s{\rm log}(d\vee n)\over n}}.$

### Proof of Theorem 3.

Theorem 2 implies

$\displaystyle\left\|\bm{\theta}_{dis}^{1}-\bm{\theta}^{*}\right\|_{2}$
$\displaystyle\lesssim\left(s\sqrt{{{\rm log}(d\vee n)\over
n}}\left\|{\bm{\theta}_{dis}^{0}}-\bm{\theta}^{*}\right\|_{2}+\sqrt{s}\left\|{\bm{\theta}_{dis}^{0}}-\bm{\theta}^{*}\right\|_{2}^{2}+\sqrt{{s{\rm
log}(d\vee n)\over N}}\right).$

This implies that after $Q\asymp{{\rm log}N\over{\rm log}n}$ iterations for
the outer loop (note if $n\geq N^{c}$ for some constant $c\in(0,1)$, $Q$ is a
constant) we obtain

$\left\|\bm{\theta}_{dis}^{Q}-\bm{\theta}^{*}\right\|_{2}\lesssim\sqrt{{s{\rm
log}(d\vee n)\over N}}.$

### Technical Lemmas

######  Lemma 1.

Given $\bm{\theta}=(\bm{\delta},\gamma)$ satisfying
$||\bm{\theta}-\bm{\theta}^{*}||_{0}\leq(2s+s_{0}+1)$ and $\gamma\geq C^{*}$,
under Assumptions 1, 2 and 3, then, with probability at least $1-(d\vee
n)^{-C}$,

$\left\|\left(\nabla^{2}\mathcal{L}_{1}\left({\bm{\theta}}\right)-\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)(\bm{\theta}-\bm{\theta}^{*})\right\|_{\infty}\lesssim\|{\bm{\theta}}-\bm{\theta}^{*}\|_{2}^{2}.$

Proof. Denote
$\bm{\theta}-\bm{\theta}^{*}=\bm{\Delta}=\left(\bm{\Delta}_{\bm{\delta}},\bm{\Delta}_{\gamma}\right)$.

  * $\bullet$

For $j=1,\dots,d+1$,

$\displaystyle{\left[\left(\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}+\bm{\Delta}\right)-\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)\bm{\Delta}\right]}_{j}$
$\displaystyle={1\over
n}\sum_{i=1}^{n}\left(1-d_{i}\right)\left(x_{ij}\mathbf{x}_{i}^{\prime}\bm{\Delta}_{\bm{\delta}}\left(h\left(\mathbf{x}_{i}^{\prime}\bm{\delta}^{*}\right)-h\left(\mathbf{x}_{i}^{\prime}\bm{\delta}\right)\right)\right.$
$\displaystyle\leq{1\over n}\sum_{i}|x_{ij}|({\bf
x}_{i}^{\prime}\mbox{\boldmath$\Delta$}_{\bm{\delta}})^{2}\lesssim\|\bm{\delta}\|_{2}^{2},$

where we use $|h^{\prime}(a)|<4.3$ for all $a\in\mathbb{R}$ (Jacobson and Zou
(2024)).

  * $\bullet$

For $j=d+2$,

$\displaystyle{\left[\left(\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}+\bm{\Delta}\right)-\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)\bm{\Delta}\right]_{j}}$
$\displaystyle={1\over
n}\sum_{i=1}^{n}d_{i}\left((\gamma)^{-2}-\left(\gamma^{*}\right)^{-2}\right)\bm{\Delta}_{\gamma}$
$\displaystyle\leq{2\over\left(C^{*}\right)^{3}}\bm{\Delta}_{\gamma}^{2},$

where the last inequality holds because of $\gamma\geq C^{*}$.

Thus, with probability of at least $1-(d\vee n)^{-C}$,

$\displaystyle\left\|\left(\nabla^{2}\mathcal{L}_{1}(\bm{\theta})-\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)\bm{\Delta}\right\|_{\infty}$
$\displaystyle\lesssim\|\bm{\Delta}\|_{2}^{2}.$

######  Lemma 2.

For ${\bm{\theta}}=({\bm{\delta}},{\gamma})$ satisfying
$||{\bm{\theta}}-\bm{\theta}^{*}||_{0}\leq(2s+s_{0}+1)$ and ${\gamma}\geq
C^{*}$, under Assumptions 1, 2, and 3, with probability at least $1-(d\vee
n)^{-C}$, we have with probability $1-(d\vee n)^{-C}$,

$\displaystyle\|\nabla{\cal L}_{1}(\mbox{\boldmath$\theta$})-\nabla{\cal
L}_{1}(\mbox{\boldmath$\theta$}^{*})-E[{\cal
L}_{1}(\mbox{\boldmath$\theta$})-\nabla{\cal
L}_{1}(\mbox{\boldmath$\theta$}^{*})]\|_{\infty}\leq C\sqrt{{s{\rm log}(d\vee
n)\over
n}}\|\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{*}\|_{2}+C\|\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{*}\|^{2}_{2}.$

Proof. We have, by Taylor’s expansion,

$\displaystyle\left(\nabla\mathcal{L}_{1}({\bm{\theta}})-\nabla\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)$
$\displaystyle=$
$\displaystyle\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\left({\bm{\theta}}-\bm{\theta}^{*}\right)+\int_{0}^{1}\left(\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}+u\left({\bm{\theta}}-\bm{\theta}^{*}\right)\right)-\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)\left({\bm{\theta}}-\bm{\theta}^{*}\right)du.$

By Lemma 1, the second term above as well as its expectation (can be shown
following exactly the same lines) is bounded by
$C\|\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{*}\|_{2}^{2}$. For the
1st term above, we first derive a bound for

$\displaystyle\left\|\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)-\nabla^{2}\mathcal{L}^{*}\left(\bm{\theta}^{*}\right)\right\|_{\max},$
(11)

where we denote $\|H\|_{\max}\doteq\max{|H_{ij}|}$ for a matrix $H$ and we
recall that by our notation
$\mathcal{L}^{*}(\bm{\theta}^{*})=E[\mathcal{L}_{1}(\bm{\theta}^{*})]$. We
analyze each entry of the empirical Hessian matrix $\nabla^{2}{\cal
L}_{1}(\mbox{\boldmath$\theta$}^{*})$.

  * $\bullet$

${\partial^{2}\over\partial\bm{\delta}_{j}\partial\bm{\delta}_{k}}\mathcal{L}_{1}(\bm{\theta})={1\over
n}\sum_{i=1}^{n}x_{ij}x_{ik}\left[d_{i}+\left(1-d_{i}\right)h\left(-\mathbf{x}_{i}^{\prime}\bm{\delta}\right)\right].$
Since $0<h(s)<1$ for all $s\in\mathbb{R}$, we see that

$|x_{ij}x_{ik}\left[d_{i}+\left(1-d_{i}\right)h\left(-\mathbf{x}_{i}^{\prime}\bm{\delta}\right)\right]|\leq|x_{ij}x_{ik}|.$

Thus,
$x_{ij}x_{ik}\left[d_{i}+\left(1-d_{i}\right)h\left(-\mathbf{x}_{i}^{\prime}\bm{\delta}\right)\right]$
is sub-exponential.

  * $\bullet$

${\partial^{2}\over\partial\bm{\delta}_{j}\partial\gamma}\mathcal{L}_{1}(\bm{\theta})={1\over
n}\sum_{i=1}^{n}-d_{i}y_{i}x_{ij}$. It is easy to see that $d_{i}y_{i}x_{ij}$
is sub-exponential.

  * $\bullet$

${\partial^{2}\over\partial\gamma^{2}}\mathcal{L}_{1}(\bm{\theta})={1\over
n}\sum_{i=1}^{n}d_{i}\left(\gamma^{-2}+y_{i}^{2}\right)$, and
$d_{i}\left(\gamma^{-2}+y_{i}^{2}\right)$ is again sub-exponential.

With the aid of the concentration inequality for sub-exponential variables,
(11) is bounded by $C\sqrt{{{\rm log}(d\vee n)\over n}}$ with probability at
least $1-(d\vee n)^{-C}$. Consequently,

$\left\|(\nabla^{2}\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)-\nabla^{2}\mathcal{L}^{*}\left(\bm{\theta}^{*}\right))({\bm{\theta}}-\bm{\theta}^{*})\right\|_{\infty}\leq
C\sqrt{{{\rm log}(d\vee n)\over n}}||{\bm{\theta}}-\bm{\theta}^{*}||_{1}\leq
C\sqrt{{s{\rm log}(d\vee n)\over n}}\|{\bm{\theta}}-\bm{\theta}^{*}\|_{2}.$

This completes the proof.

######  Lemma 3.

Given $\bar{\bm{\theta}}=(\bar{\bm{\delta}},\bar{\gamma})$ satisfying
$||\bar{\bm{\theta}}-\bm{\theta}^{*}||_{0}\leq(2s+s_{0}+1)$ and
$\bar{\gamma}\geq C^{*}$, under Assumptions 1, 2, and 3, with probability at
least $1-(d\vee n)^{-C}$, we have,

$\displaystyle\left\|\nabla\widetilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right\|_{\infty}$
$\displaystyle\lesssim\sqrt{{s{\rm log}(d\vee n)\over
n}}\left\|\overline{\bm{\theta}}-\bm{\theta}^{*}\right\|_{2}+\|\overline{\bm{\theta}}-\bm{\theta}\|_{2}^{2}+\sqrt{{{\rm
log}(d\vee n)\over N}}.$

Proof.

$\displaystyle\nabla\tilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)=$
$\displaystyle\nabla\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)-\nabla\mathcal{L}_{1}(\bar{\bm{\theta}})+\mathcal{L}_{N}(\bar{\bm{\theta}})$
$\displaystyle=$
$\displaystyle\left(\nabla\mathcal{L}_{N}(\bar{\bm{\theta}})-\nabla\mathcal{L}_{N}\left(\bm{\theta}^{*}\right)\right)-\left(\nabla\mathcal{L}_{1}(\bar{\bm{\theta}})-\nabla\mathcal{L}_{1}\left(\bm{\theta}^{*}\right)\right)+\nabla\mathcal{L}_{N}\left(\bm{\theta}^{*}\right).$

By Lemma 2, with probability at least $1-(d\vee n)^{-C}$, we have

$\displaystyle\left\|\nabla\widetilde{\mathcal{L}}\left(\bm{\theta}^{*}\right)\right\|_{\infty}\lesssim$
$\displaystyle\sqrt{{s{\rm log}(d\vee n)\over
n}}\left\|\overline{\bm{\theta}}-\bm{\theta}^{*}\right\|_{2}+\|\overline{\bm{\theta}}-\bm{\theta}\|_{2}^{2}+\left\|\nabla\mathcal{L}_{N}\left(\bm{\theta}^{*}\right)\right\|_{\infty}.$

Recall that $0<|g^{\prime}(s)|<1$, and thus
$|g(-\mathbf{x}^{\prime}\bm{\delta})-g(0)|\leq|{\bf
x}^{\prime}\mbox{\boldmath$\delta$}|$, implying
$g(-\mathbf{x}^{\prime}\bm{\delta})$ is sub-Gaussian. Thus using the same
arguments as in the proof of Lemma 2, it is easy to see, based on the
expression (6) for the gradient, that the components of
$\mathcal{L}_{N}(\bm{\theta}^{*})$ are all sub-exponential. Thus
$\left\|\nabla\mathcal{L}_{N}\left(\bm{\theta}^{*}\right)\right\|_{\infty}\lesssim\sqrt{{{\rm
log}(d\vee n)\over N}}$ holds with probability at least $1-(d\vee n)^{-C}$.

## References

  * (1)
  * Belloni and Chernozhukov (2011) Belloni, A., and Chernozhukov, V. (2011), “l1-penalized quantile regression in high-dimensional sparse models,” The Annals of Statistics, 39(1), 82–130.
  * Blumensath and Davies (2009) Blumensath, T., and Davies, M. E. (2009), “Iterative hard thresholding for compressed sensing,” Applied and Computational Harmonic Analysis, 27(3), 265–274.
  * Buckley and James (1979) Buckley, J., and James, I. (1979), “Linear regression with censored data,” Biometrika, 66(3), 429–436.
  * Fan and Li (2001) Fan, J., and Li, R. (2001), “Variable selection via nonconcave penalized likelihood and its oracle properties,” Journal of the American Statistical Association, 96(456), 1348–1360.
  * Foucart (2011) Foucart, S. (2011), “Hard thresholding pursuit: an algorithm for compressive sensing,” SIAM Journal on Numerical Analysis, 49(6), 2543–2563.
  * Garg and Khandekar (2009) Garg, R., and Khandekar, R. (2009), Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property,, in Proceedings of the 26th annual international conference on machine learning, pp. 337–344.
  * Jacobson and Zou (2024) Jacobson, T., and Zou, H. (2024), “High-dimensional censored regression via the penalized Tobit likelihood,” Journal of Business & Economic Statistics, 42(1), 286–297.
  * Jain et al. (2014) Jain, P., Tewari, A., and Kar, P. (2014), On iterative hard thresholding methods for high-dimensional M-estimation,, in Advances in Neural Information Processing Systems.
  * Ji et al. (2023) Ji, Y., Scutari, G., Sun, Y., and Honnappa, H. (2023), “Distributed sparse regression via penalization,” Journal of Machine Learning Research, 24(272), 1–62.
  * Jordan et al. (2019) Jordan, M. I., Lee, J. D., and Yang, Y. (2019), “Communication-efficient distributed statistical inference,” Journal of the American Statistical Association, 114(526), 668–681.
  * Lee et al. (2017) Lee, J. D., Liu, Q., Sun, Y., and Taylor, J. E. (2017), “Communication-efficient sparse regression,” Journal of Machine Learning Research, 18(5), 1–30.
  * Müller and van de Geer (2016) Müller, P., and van de Geer, S. (2016), “Censored linear model in high dimensions,” TEST, 25, 75–92.
  * Nedic and Ozdaglar (2009) Nedic, A., and Ozdaglar, A. (2009), “Distributed subgradient methods for multi-agent optimization,” IEEE Transactions on Automatic Control, 54(1), 48–61.
  * Olsen (1978) Olsen, R. J. (1978), “Note on the uniqueness of the maximum likelihood estimator for the Tobit model,” Econometrica, 46(5), 1211–1215.
  * Sampford (1953) Sampford, M. R. (1953), “Some inequalities on Mill’s ratio and related functions,” The Annals of Mathematical Statistics, 24(1), 130 – 132.
  * Soret et al. (2018) Soret, P., Avalos, M., Wittkop, L., Commenges, D., and Thiébaut, R. (2018), “Lasso regularization for left-censored Gaussian outcome and high-dimensional predictors,” BMC Medical Research Methodology, 18(1), 159.
  * Tobin (1958) Tobin, J. (1958), “Estimation of relationships for limited dependent variables,” Econometrica, 26(1), 24–36.
  * Wang et al. (2023) Wang, Y., Lu, W., and Lian, H. (2023), “Best subset selection for high-dimensional non-smooth models using iterative hard thresholding,” Information Sciences, 625, 36–48.
  * Zhang et al. (2013) Zhang, Y., Duchi, J. C., and Wainwright, M. J. (2013), “Communication-efficient algorithms for statistical optimization,” Journal of Machine Learning Research, 14(68), 3321–3363.
  * Zhou and Liu (2016) Zhou, X., and Liu, G. (2016), “LAD-Lasso variable selection for doubly censored median regression models,” Communications in Statistics - Theory and Methods, 45(12), 3658–3667.


1

  * June 2020

# Modelling lattices of organic polaritons

Thomas J Sturges1 and Magdalena Stobińska1 1Institute of Theoretical Physics,
University of Warsaw, ul. Pasteura 5, 02-093, Warsaw, Poland
<EMAIL_ADDRESS>1Institute of Theoretical Physics, University of Warsaw,
ul. Pasteura 5, 02-093, Warsaw, Poland<EMAIL_ADDRESS>
###### Abstract

Microcavity-polaritons in two-dimensional lattice geometries have been used to
study a wide range of interesting physics. Meanwhile, organic materials have
shown great promise on the road towards polaritonic devices, as the strong
binding energy of their Frenkel excitons permits room temperature condensation
and lasing. Whilst there are theoretical treatments of the condensation
processes in planar organic microcavities, models of lattice geometries are
lacking. Here, we introduce a model for describing the dynamics of lattices of
zero-dimensional organic microcavities, where the dominant condensation
mechanism involves the emission of a vibrational phonon. The vibronic
resonance provides a tool for targeted condensation in a particular eigenmode
of the system, which we highlight by examining a dimer and a dimerised chain.
For the dimer, we observe a double resonance in the condensation efficiency
that arises from tuning the condensate-reservoir detuning in to resonance with
either the symmetric or antisymmetric mode. This selective condensation was
also exploited in the dimerised chain, to condense the system into the
topological edge states under homogeneous pumping of all cavities. We also
showed an interesting signature of nonreciprocal transport when pumping a
single cavity in the chain, where the direction of propagation depends on the
sublattice being pumped.

Keywords: polaritons, organic, lattice, transport, excitons, vibrons

##  1 Introduction

Microcavity-polaritons hold much promise for the next generation of
optoelectronic devices [1]; including low-threshold lasers and light-emitters
[2], optimisation machines [3], optical circuitry and logic gates [4], as well
as a platform for quantum simulations and quantum computing [5]. They are
quasiparticles that arise from the strong light-matter coupling between
excitons in quantum wells or wires, and confined photonic modes in optical
microcavities [6, 7]. They can form polariton condensates [8] — a coherent
state in which a single mode is macroscopically occupied — in analogy to
atomic Bose-Einstein condensates [9]. In addition, modern fabrication
techniques permit fine control over the lateral confinement and potential
landscape of polaritons. This has enabled the exploration of fascinating
phenomena in lattices of polaritons [10], such as exciton-polariton
topological insulators [11], Dirac cones in polaritonic graphene [12], and
Josephson junctions [13].

Much of the pioneering work in polaritonics employed inorganic materials and
cryogenic temperatures [8]; this is because the temperature must be lower than
the binding energy of the Wannier excitons in these systems. Organic materials
are an attractive alternative as they allow room temperature condensation [14]
due to the much larger binding energy of their Frenkel excitons [15]; in
addition their processing and fabrication can also be much simpler [16].
Organic materials feature high-energy vibrational modes, which enables the
efficient condensation of polaritons when matching a vibrational resonance
[17, 18]. This opens the interesting perspective for selectively condensing
into different eigenmodes of polariton lattices [19]. Whilst there are
theoretical treatments of vibron-mediated condensation in planar microcavities
[20, 21, 22, 23, 24], a model that treats this phenomenon in lattices of
polaritons is lacking.

In this article, we introduce a general model for simulating organic exciton-
polariton condensates in lattices. Specifically we consider lattices of zero-
dimensional microcavities, each of which hosts a localised condensate mode,
that hybridise due to evanescent coupling. A nonresonant laser creates a
reservoir of excitons that can relax to the condensate by emitting a vibron;
this is a very efficient process so long as the vibronic energy is similar to
the energy difference between the reservoir and a condensate mode. This
introduces an energy-selective mechanism, in contrast to nonresonant
excitation in inorganic systems. This can be exploited to selectively populate
one or more eigenstates of the lattice and should allow the observation of
novel dynamics. We highlight the key features of this new model with a few
examples. After demonstrating the basics with a single zero-dimensional
cavity, we treat the simplest coupled system, a dimer, and observe how either
the symmetric or anti-symmetric mode can be selectively populated by
modulating the detuning. Then, we examine the dimerised chain and show that a
homogeneous pumping of the chain can result in exciting the topological mid
gap states, or alternatively the bulk states, by controlling the detuning.
Finally, we show a signature of nonreciprocal transport that emerges when
driving just one site in the middle of the chain. Taken together, these
results exemplify some of the new interesting physics to be found in lattices
of organic polaritons that feature different dynamics to their inorganic
counterparts.

The rest of this article is organised as follows. In section 2 we introduce
the analytical model; we describe the Hamiltonian and derive the mean-field
equations. Then in section 3 we demonstrate the interesting features of these
equations for simulating lattices of organic polaritons, by exploring a few
examples; a single cavity, a dimer, and a dimerised chain where we show the
formation of edge states and a signature of nonreciprocal transport. We
discuss the results in section 4, and explain more details of our methods in
section 5.

##  2 Analytical model

Zero-dimensional microcavities can host localised polariton condensates. For
example, one can carve a defect microcavity into a planar Distributed-Bragg-
Reflector by Focused-Ion Beam milling; this leads to a discretisation of the
lower polariton branch into localised modes within the defect microcavity,
with the ground state well separated from higher energy states and featuring a
flat dispersion [18]. For such systems it is reasonable for us to effectively
describe the condensate of the microcavity as a single bosonic mode. This
allows us to treat the dynamics of coupled arrays of such cavities with a
lattice model. Thus we introduce a Hamiltonian for a general array as

$H=\sum_{m}\left(\omega_{\mathrm{C}}\psi_{m}^{\dagger}\psi_{m}+\omega_{\mathrm{R}}a_{m}^{\dagger}a_{m}+\omega_{\mathrm{V}}b_{m}^{\dagger}b_{m}\right)+H_{\mathrm{I}}+H_{\mathrm{F}},$
(1)

where the energy and corresponding creation operator for the condensate are
denoted by $\omega_{\mathrm{C}}$ and $\psi_{m}^{\dagger}$, for the excitonic
reservoir by $\omega_{\mathrm{R}}$ and $a_{m}^{\dagger}$, and for the vibrons
by $\omega_{\mathrm{V}}$ and $b_{m}$; in the case of the $m$th microcavity. By
fabricating two (or more) microcavities with a small overlap, there will be
evanescent coupling between the condensates, which can be described by

$H_{\mathrm{I}}=\frac{1}{2}\sum_{mn}\Omega_{mn}\left(\psi_{m}^{\dagger}\psi_{n}+\psi_{n}^{\dagger}\psi_{m}\right),$
(2)

where $\Omega_{mn}=\Omega_{nm}$ is the coupling energy between the $m$th and
$n$th cavity. In organic microcavities, scattering from the excitonic
reservoir to the polaritonic states are often facilitated by an electronic
transition that is not strongly coupled [16, 20, 17, 18]. Here we consider a
vibron-mediated relaxation which can be modelled with a Frölich-type
Hamiltonian

$H_{\mathrm{F}}=g\sum_{m}\left(a_{m}b_{m}^{\dagger}\psi_{m}^{\dagger}+a_{m}^{\dagger}b_{m}\psi_{m}\right),$
(3)

where $g$ is the interaction constant between excitons and vibrons. We
describe particle decay and a nonresonant pump via a master equation (see
Methods), from which we can derive mean-field equations for the number of
reservoir excitons $N_{m}=\langle a_{m}^{\dagger}a_{m}\rangle$, the energy
flow $J_{mn}=\langle a_{m}b_{m}^{\dagger}\psi_{n}^{\dagger}\rangle$, the flow
of condensate polaritons $\rho_{mn}=\langle\psi_{m}^{\dagger}\psi_{n}\rangle$,
as well as the number of condensate polaritons $\rho_{m}=\rho_{mm}$. The
equations read

$\displaystyle\hbar\frac{\text{d}\rho_{mn}}{\text{d}t}$
$\displaystyle=-\gamma_{\mathrm{C}}\rho_{mn}+\text{i}g\left(J_{mn}^{*}-J_{nm}\right)+\text{i}\sum_{\ell=1}^{M}\left(\Omega_{m\ell}\rho_{n\ell}^{*}-\Omega_{n\ell}\rho_{m\ell}\right),$
$\displaystyle\hbar\frac{\text{d}N_{m}}{\text{d}t}$
$\displaystyle=-\gamma_{R}N_{m}+P_{m}-\text{i}g\left(J_{mm}^{*}-J_{mm}\right),$
(4) $\displaystyle\hbar\frac{\text{d}J_{mn}}{\text{d}t}$
$\displaystyle=-\left(\text{i}\delta+\frac{\Gamma}{2}\right)J_{mn}+\text{i}gN_{m}^{R}\left(\delta_{mn}+\rho_{nm}\right)+\text{i}\sum_{\ell=1}^{M}\Omega_{n\ell}J_{m\ell},$

where $\gamma_{\mathrm{C}}$, $\gamma_{\mathrm{R}}$, $\gamma_{\mathrm{V}}$, and
$\Gamma=\gamma_{\mathrm{C}}+\gamma_{\mathrm{R}}+\gamma_{\mathrm{V}}$ are the
inverse decay rates of the condensate, reservoir, vibrons, and their sum
respectively; $P_{m}$ is the pumping rate to the reservoir of the $m$th
cavity, and
$\delta=\omega_{\mathrm{R}}-(\omega_{\mathrm{C}}+\omega_{\mathrm{V}})$ is the
detuning between the reservoir and the sum of the condensate and vibron
energies. Thus we have a system of $M(2M+1)$ equations and variables (where
$M$ is the number of lattice sites) which we can solve for the dynamics. We
now gain some insight and intuition with these equations by exploring some
simple models, and also observe a few interesting results.

##  3 Results

###  3.1 Single Cavity

In this section only, we assume that the energy flows can be adiabatically
eliminated; in which case equation 4 simplifies to

$\displaystyle\hbar\frac{\text{d}N}{\text{d}t}$
$\displaystyle=-\gamma_{\mathrm{R}}N+P(t)-GN\left(\rho+1\right),$ (5)
$\displaystyle\hbar\frac{\text{d}\rho}{\text{d}t}$
$\displaystyle=-\gamma_{\mathrm{C}}\rho+GN\left(\rho+1\right),$

where

$G=\frac{g^{2}\Gamma}{\delta^{2}+(\Gamma/2)^{2}},$ (6)

is the condensation efficiency. These equations clearly show how the
efficiency with which reservoir excitons transfer into the condensate depends
on the vibronic transition being in resonance with the energy difference
between the reservoir and condensate; i.e. $G$ is a maximum for $\delta=0$.
Figure 1 shows a close correspondence between the value of $G$ and the number
of polaritons, obtained from a steady-state solution of equation 5 (with
$P(t)=P_{0}$). We note that the condensation process here is quite distinct to
both resonant and non-resonant laser excitation of inorganic systems. We have
a condensation mechanism which is incoherent and phase non-preserving (like
non-resonant pumping of inorganic materials); but which, due to the vibronic
resonance, can be used to selectively condense at a specific energy (like
resonant pumping of inorganic materials).

Figure 1: Single Cavity. The number of condensate polaritons $\rho$, and the
parameter $G$ that specifies the efficiency of vibron-mediated condensation.
The parameters are given in Methods.

###  3.2 Dimer

The simplest example of an interacting array of condensates is a dimer; which
can also be thought of as the building block (the unit cell) of more
complicated bipartite lattices. The hybridisation between the two condensates
leads to a symmetric and anti-symmetric mode with an energy splitting equal to
the coupling constant. By solving equation 4 for the steady-state solutions
(where $M=2$ and $\Omega_{12}\equiv\Omega_{0}$) we obtain solutions for the
number of condensate polaritons as a function of detuning; see Figure 2. When
pumping both sites equally (Figure 2(a)) the population on both sites is,
unsurprisingly, the same. Moreover, we see a double resonance corresponding to
peaks at the symmetric and anti-symmetric modes. When pumping just the left
site (Figure 2(b)) we see a curious effect. For large detunings, far away from
resonance, the pumped left (red) site has a larger population than the non-
pumped right (black) site. However for detunings in-between the two
(symmetric/anti-symmetric) resonances the population flows from the pumped to
the non-pumped site and builds up there. We now look at a one-dimensional
chain constructed from this building block, the dimerised chain.

(a) Both sites pumped the same

(b) Only left ($m=1$) is pumped

Figure 2: Double cavity. The number of condensate polaritons in the left
cavity ($m=1$, red line) and right cavity ($m=2$, black dashed line) as a
function of reservoir-condensate detuning $\delta$. (a) Both sites are pumped
the same; $P_{1}(t)=P_{2}(t)=P_{0}$. (b) Only the left cavity is pumped;
$P_{1}(t)=P_{0}$, $P_{2}(t)=0$. The parameters are given in Methods.

###  3.3 The topological midgap states of a dimerised chain

We now consider a one-dimensional dimerised chain, with alternating coupling
strengths between nearest neighbours. As such the condensate interaction
Hamiltonian (equation 2) embodies the Su-Schrieffer-Heeger model which
features topologically protected midgap states that are exponentially
localised at the edges. We encode the geometry by setting the nearest
neighbour couplings in equation 4 to $\Omega_{m,m\pm
1}=\Omega_{0}\pm(-1)^{m}\delta_{\Omega}$ and all other couplings to zero. We
consider a homogeneous pumping profile where all sites are pumped equally.

(a) Eigenvalues of $\mathcal{H}$

(b) Eigenvector of the midgap state

Figure 3: Eigenstates of the dimerised chain. (a) The eigenvalues of the
condensate Hamiltonian $\mathcal{H}$, with the fill colour of the circles
coded by the corresponding value of $G$ (equation 6) and (b) the eigenvector
corresponding to one of the midgap states (Eigenvalue no. 10). The parameters
are given in Methods.

Before simulating the dynamics, let us gain some insight into the system by
finding the eigenstates of the condensate interaction Hamiltonian. We can
write equation 2 in matrix form as
$H_{\mathrm{I}}=\boldsymbol{\hat{\psi}}^{\dagger}\mathcal{H}\boldsymbol{\hat{\psi}}$,
where $\boldsymbol{\hat{\psi}}=(\psi_{1},\psi_{2},\dots,\psi_{M})^{T}$; we
show the eigenvalues of $\mathcal{H}$ in Figure 3(a), and the eigenvector of
one of the midgap states in Figure 3(b). In addition, the fill colour of the
circles that mark the eigenvalues denotes the single cavity condensation
efficiency $G$, equation 6. In this way, we can see that we have chosen
parameters such that the vibronic transition is in resonance with the midgap
states (which also happens to be the energy of a single isolated cavity). As
such the condensation process should predominantly excite the midgap states,
and we expect to see the formation of edge states similar to the pure
eigenvectors; see Figure 3(b).

To simulate the dynamics in the time domain we use the Runga-Kutta method of
fourth order (see Methods). We start the simulation with an empty system,
$\rho_{mn}=N_{m}=J_{mn}=0$, and switch on a constant homogeneous pump
$P_{m}(t)=P_{0}$. As we can see in Figure 4, the system condenses to achieve a
steady state solution that features a strong signature of the midgap edge
states. Of course the condensation efficiency $G$ is not a delta-function, but
features a broadening due to the total decay rate $\Gamma$; and thus the
steady-state solution also features a weaker mix of some of the other bulk
states, which acts to wash out the pure edge state. Nonetheless a pronounced
edge state is indeed witnessed.

In order to test the robustness of this result, we performed stochastic
simulations of the same system with a realistic disorder in both the value of
$\omega_{\mathrm{C}}$ at each site and the coupling strengths $\Omega_{mn}$.
In both cases the disorder was sampled from a normal distribution with
standard deviations $\sigma_{\mathrm{C}}=0.04\Omega_{0}$ and
$\sigma_{\Omega}=0.01\Omega_{0}$ respectively; see Methods for details. As can
be seen from the error bars in Figure 4(b), the key features of the edge
states are not washed out by disorder.

(a) Simulation of equation 4

(b) Bar chart of final time step

Figure 4: Condensation in the edge states. (a) A time-domain simulation of
equation 4 for the dimerised chain, and (b) a bar chart of the final time
step. The error bars denote the standard deviation of the corresponding
stochastic simulations. The parameters are given in Methods.

###  3.4 Nonreciprocal transport

In the previous section we showed how the system selectively excites the edge
states under a homogeneous pumping of all cavities. Of course we are free to
choose the pumping pattern as we like; and so we can explore what interesting
effects might arise in other circumstances. As one example, we investigate
what happens when pumping just one cavity in the middle of a chain with $M=40$
sites. First, we set $P_{m}(t)=P_{0}\delta_{m,20}$ and simulate the dynamics
of equation 4. As shown in Figure 5(a), the polaritons propagate predominantly
to one side of the excitation spot. If we now move the pump to a neighbouring
chain, such that $P_{m}(t)=P_{0}\delta_{m,21}$, the polaritons now propagate
in the opposite direction; see Figure 5(b). As such the system is displaying a
type of nonreciprocal transport, where the propagation direction is determined
by which sublattice is being pumped.

(a) Pumping site no. 5

(b) Pumping site no. 6

Figure 5: Nonreciprocal transport. The parameters are given in Methods.

##  4 Discussion

We have introduced a model for describing the dynamics of organic polaritons
trapped in lattices of zero-dimensional microcavities; which interact due to
evanescent coupling, and condense by emitting a vibrational phonon. Our
treatment assumes that condensation occurs in one mode per microcavity, but
this can be easily generalised to multimode cavities. Such a lattice model
cannot capture physics related to the details of specific trapping geometries,
such as the spatial density patterns that arise when condensing in different
order orbital modes. Nonetheless, it is an important step in the theoretical
treatment of condensation processes in lattice of organic polaritons.

Our results demonstrate the need for methods in modelling lattices of organic
polaritons, other than those typically used in inorganic polaritons, such as
generalised complex Gross-Pitaevskii equations which do not capture the new
physics that arises due to the vibron-mediated condensation. As such, our
model contributes to the theoretical toolbox, and should be readily applied to
other organic lattices.

##  5 Methods

The mean-field equations 4 are obtained from the Hamiltonian in equation 1 and
the master equation

$\displaystyle\hbar\frac{\text{d}\hat{\boldsymbol{\rho}}}{\text{d}t}=\text{i}\left[\hat{\boldsymbol{\rho}},H\right]$
$\displaystyle+\sum_{n}\left[\gamma_{0}\mathcal{L}(\psi_{n})+\gamma_{R}\mathcal{L}(a_{n})+P_{n}(t)\left(\mathcal{L}(a_{n})+\mathcal{L}(a_{n}^{\dagger})\right)\right]$
(7)
$\displaystyle+\sum_{n}\gamma_{b}\left[(1+n_{\text{th}})\mathcal{L}(b_{n})+n_{\text{th}}\mathcal{L}(b_{n}^{\dagger})\right],$

where $\hat{\boldsymbol{\rho}}$ is the density matrix, $n_{\mathrm{th}}$ is
the thermal distribution of vibrons, and
$\mathcal{L}(a)=(1/2)(2a\hat{\boldsymbol{\rho}}a^{\dagger}-a^{\dagger}a\hat{\boldsymbol{\rho}}-\hat{\boldsymbol{\rho}}a^{\dagger}a)$.
We use the identity
$\frac{\mathrm{d}\langle{\hat{A}}\rangle}{\mathrm{d}t}=\mathrm{Tr}(\frac{\mathrm{d}\hat{\boldsymbol{\rho}}}{\mathrm{d}t}\hat{A})$
which leads to the mean-field equations. We make two approximations; we assume
that at any time the number, and thus flow, of vibrons is negligible $\langle
b_{m}^{\dagger}b_{n}\rangle\ll 1$, and we discard nonresonant processes. With
these assumptions we obtain the mean-field equations 4.

| $\omega_{\mathrm{C}}$ | $\omega_{\mathrm{R}}$ | $\omega_{\mathrm{V}}$ | $\gamma_{\mathrm{C}}$ | $\gamma_{\mathrm{R}}$ | $\gamma_{\mathrm{V}}$ | $\Omega_{0}$ | $\delta_{\Omega}$ | $P_{0}$ | $g$ | $\sigma_{\mathrm{C}}$ | $\sigma_{\Omega}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Figure 1 | 2526 | $\omega_{\mathrm{R}}$ | 194 | 4 | 0.7 | 2.9 | - | - | 30 | 0.5 | - | -
Figure 2 | 2526 | $\omega_{\mathrm{R}}$ | 194 | 4 | 0.7 | 2.9 | 7 | - | 30 | 0.5 | - | -
Figure 3 | 2526 | 2720 | 194 | 4 | 0.7 | 2.9 | 33 | 7 | - | - | - | -
Figure 4 | 2526 | 2720 | 194 | 4 | 0.7 | 2.9 | 33 | 7 | 30 | 0.5 | 1.32 | 3.3
Figure 5 | 2526 | 2720 | 194 | 4 | 0.7 | 2.9 | 33 | 7 | 30 | 0.5 | - | -
Table 1: Parameters used. In this table we list the values of the parameters
used in all of the Figures in meV.

To obtain the steady-state solutions we set the time derivatives to zero in
equation 5 for the single cavity, and in equation 4 with $M=2$ for the dimer.
We can then solve the resulting equations analytically in the case of a single
cavity, and numerically with Mathematica [25] for the dimer. To simulate the
dynamics of equation 4 we use XMDS2 [26] which is an open source package for
numerically integrating initial value problems, including differential
equations. We use the solver which employs the Runga-Kutta method of fourth
order. To perform the stochastic simulations of Figure 4(b) we modify equation
4 with the substitutions
$\omega_{\mathrm{C}}\xrightarrow[]{}\omega_{\mathrm{C}}+\mathcal{R}_{m}(\sigma_{\mathrm{C}})$
and $\Omega_{mn}\xrightarrow[]{}\Omega_{mn}+\mathcal{R}_{mn}(\sigma_{\Omega})$
where $\mathcal{R}_{i}(\sigma)$ is a random number sampled from a normal
distribution with variance $\sigma^{2}$. We then perform 100 simulations of
different instances of the disorder; the bar chart in Figure 4(b) show the
average, and the error bars show the standard deviation of simulation results.
The source code is available for download; see the Supplementary Material. In
Table 1 we gather the parameters used in all figures in the manuscript.

## References

  * [1] Daniele Sanvitto and Stéphane Kéna-Cohen  “The road towards polaritonic devices”  In _Nature Materials_ 15.10, 2016, pp. 1061–1073  URL: https://doi.org/10.1038/nmat4668
  * [2] Ayan Das et al.  “Room Temperature Ultralow Threshold GaN Nanowire Polariton Laser”  In _Phys. Rev. Lett._ 107 American Physical Society, 2011, pp. 066405  DOI: 10.1103/PhysRevLett.107.066405
  * [3] Natalia G. Berloff et al.  “Realizing the classical XY Hamiltonian in polariton simulators”  In _Nature Materials_ 16.11, 2017, pp. 1120–1126  URL: https://doi.org/10.1038/nmat4971
  * [4] T… Liew, A.. Kavokin and I.. Shelykh  “Optical Circuits Based on Polariton Neurons in Semiconductor Microcavities”  In _Phys. Rev. Lett._ 101 American Physical Society, 2008, pp. 016402  DOI: 10.1103/PhysRevLett.101.016402
  * [5] Dimitris G Angelakis  “Quantum Simulations with Photons and Polaritons”  Springer International Publishing, 2017
  * [6] J.. Hopfield  “Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals”  In _Phys. Rev._ 112 American Physical Society, 1958, pp. 1555–1567  DOI: 10.1103/PhysRev.112.1555
  * [7] Guillaume Malpuech Alexey V. and Fabrice P. Laussy  “Microcavities”  Oxford Univeristy Press, 2017
  * [8] Tim Byrnes, Na Young Kim and Yoshihisa Yamamoto  “Exciton-polariton condensates”  In _Nature Physics_ 10.11, 2014, pp. 803–813  URL: https://doi.org/10.1038/nphys3143
  * [9] M.. Anderson et al.  “Observation of Bose-Einstein Condensation in a Dilute Atomic Vapor”  In _Science_ 269.5221, 1995, pp. 198–201  DOI: 10.1126/science.269.5221.198
  * [10] C Schneider et al.  “Exciton-polariton trapping and potential landscape engineering”  In _Reports on Progress in Physics_ 80.1 IOP Publishing, 2016, pp. 016503  DOI: 10.1088/0034-4885/80/1/016503
  * [11] S. Klembt et al.  “Exciton-polariton topological insulator”  In _Nature_ 562.7728, 2018, pp. 552–556  URL: https://doi.org/10.1038/s41586-018-0601-5
  * [12] T. Jacqmin et al.  “Direct Observation of Dirac Cones and a Flatband in a Honeycomb Lattice for Polaritons”  In _Phys. Rev. Lett._ 112 American Physical Society, 2014, pp. 116402  DOI: 10.1103/PhysRevLett.112.116402
  * [13] K.. Lagoudakis et al.  “Coherent Oscillations in an Exciton-Polariton Josephson Junction”  In _Phys. Rev. Lett._ 105 American Physical Society, 2010, pp. 120403  DOI: 10.1103/PhysRevLett.105.120403
  * [14] S. Ké’na-Cohen and S.. Forrest  “Room-temperature polariton lasing in an organic single-crystal microcavity”  In _Nature Photonics_ 4.6, 2010, pp. 371–375  URL: https://doi.org/10.1038/nphoton.2010.86
  * [15] V.. G. Bassani  “Electronic Excitations in Organic Based Nanostructures, Volume 31”, Thin Films and Nanostructures 31  Academic Press, 2003  URL: http://gen.lib.rus.ec/book/index.php?md5=5301F7D895895B3C0E11973FEA461977
  * [16] Jonathan Keeling and Stéphane Kéna-Cohen  “Bose-Einstein Condensation of Exciton-Polaritons in Organic Microcavities”  In _Annual Review of Physical Chemistry_ 71.1, 2020, pp. 435–459  DOI: 10.1146/annurev-physchem-010920-102509
  * [17] Johannes D. Plumhof et al.  “Room-temperature Bose-Einstein condensation of cavity exciton-polaritons in a polymer”  In _Nature Materials_ 13.3, 2014, pp. 247–252  URL: https://doi.org/10.1038/nmat3825
  * [18] Fabio Scafirimuto et al.  “Room-Temperature Exciton-Polariton Condensation in a Tunable Zero-Dimensional Microcavity”  In _ACS Photonics_ 5.1 American Chemical Society, 2018, pp. 85–89  DOI: 10.1021/acsphotonics.7b00557
  * [19] Fabio Scafirimuto et al.  “Tunable exciton-polariton condensation in a two-dimensional Lieb lattice at room temperature”  In _Communications Physics_ 4.1, 2021, pp. 39  URL: https://doi.org/10.1038/s42005-021-00548-w
  * [20] L. Mazza, L. Fontanesi and G.. La Rocca  “Organic-based microcavities with vibronic progressions: Photoluminescence”  In _Phys. Rev. B_ 80 American Physical Society, 2009, pp. 235314  DOI: 10.1103/PhysRevB.80.235314
  * [21] L. Mazza, S. Kéna-Cohen, P. Michetti and G.. La Rocca  “Microscopic theory of polariton lasing via vibronically assisted scattering”  In _Phys. Rev. B_ 88 American Physical Society, 2013, pp. 075321  DOI: 10.1103/PhysRevB.88.075321
  * [22] Artem Strashko, Peter Kirton and Jonathan Keeling  “Organic Polariton Lasing and the Weak to Strong Coupling Crossover”  In _Phys. Rev. Lett._ 121 American Physical Society, 2018, pp. 193601  DOI: 10.1103/PhysRevLett.121.193601
  * [23] Luis A. Martínez-Martínez, Elad Eizner, St[é]phane Kéna-Cohen and Joel Yuen-Zhou  “Triplet harvesting in the polaritonic regime: A variational polaron approach”  In _The Journal of Chemical Physics_ 151.5, 2019, pp. 054106  DOI: 10.1063/1.5100192
  * [24] Justyna A. Ćwik, Sahinur Reja, Peter B. Littlewood and Jonathan Keeling  “Polariton condensation with saturable molecules dressed by vibrational modes”  In _EPL (Europhysics Letters)_ 105.4 IOP Publishing, 2014, pp. 47009  DOI: 10.1209/0295-5075/105/47009
  * [25] Wolfram Research, Inc.  “Mathematica, Version 12.0”  URL: https://www.wolfram.com/mathematica
  * [26] Graham R. Dennis, Joseph J. Hope and Mattias T. Johnsson  “XMDS2: Fast, scalable simulation of coupled stochastic partial differential equations”  In _Computer Physics Communications_ 184.1, 2013, pp. 201–208  DOI: https://doi.org/10.1016/j.cpc.2012.08.016
  * [27] , https://github.com/tomsturges/organic-polariton-lattices
  * [28] , https://github.com/jfeist/QuantumAlgebra.jl

T.S. and M.S. were supported by the Foundation for Polish Science “First Team”
project No. POIR.04.04.00-00-220E/16-00 (originally: FIRST TEAM/2016-2/17).
M.S. was also supported by the National Science Centre ‘Sonata Bis’ project
No. 2019/34/E/ST2/00273. We are very grateful to Rainer F. Mahrt, Thilo
Stöferle, and Darius Urbonas from IBM Research Zurich for many helpful
discussions throughout this project. We would also like to thank Anton
Zasedatelev (Skoltech) and Evgeny S. Andrianov (MIPT) for insightful
discussions at the beginning of this project. Supplementary Material:
Modelling lattices of organic polaritons Thomas J Sturges1 and Magdalena
Stobińska1 1

  * June 2020

##  1 More information on the code used for our simulations

The source code used for the simulations in the main text can be downloaded
from GitHub [27]. This custom script (org-pol-lat.xmds) was written using the
software package XMDS2 [26]. As explained in the Methods section, we use the
Runga-Kutta method of fourth order to solve equation (4) in the main text.
Here we provide some additional guidance on the source code. We encourage
anybody who is interested in adapting the code to their own needs to get in
touch with us.

Table 2 provides a translation between the LaTeX symbols used in the main text
and the plain text variables used in the XMDS2 language. Let us now run
through the main parts of the source code (please also consult the XMDS2
documentation). We use units of $meV$ and $ps$. We define our parameters
inside <arguments>, and inside <globals> we have the helpful function Kron(m,
n) to represent the Kronecker delta function $\delta_{mn}$. We represent the
subscripts $m$ and $n$ as dimensions in our simulation (so-called
<transverse_dimensions>), in addition to the main temporal dimension (the
<propagation_dimension>). The number of lattice sites is specified by M.

LaTeX | XMDS2
---|---
$\omega_{\mathrm{X}}$, $\gamma_{\mathrm{X}}$, $g$, $\Gamma$, $\Omega_{0}$, $P_{0}$, $\hbar$ | omegaX, gammaX, g, Gamma, Omega0, P0, hbar
$\rho_{mn}$, $J_{mn}$, $\Omega_{mn}$ | rho(m => m, n => n), J(…), Omega(…)
$N_{m}$, $P_{m}$, $\delta$ | NR(m => m), P(…), delta(…)
Table 2: Translation between the LaTeX symbols of the analytical model and the
plain text variables used in the code.

The variables we are solving for ($\rho_{mn}$, $J_{mn}$ and $N_{m}$) are
specified inside the <vector> named ‘variables’, which are indexed as e.g.
$J_{xy}\xrightarrow[]{}$J(m => x, n => y). We use <noise_vector> to provide
random numbers that introduce disorder for the the couplings $\Omega$ and the
detunings $\delta$, which are specified in the vectors ‘OmegaVec’ and
‘deltaVec’. The strength of the disorder is specified by sigmaOmega and
sigmadelta respectively.

The differential equations are specified inside <integrate> where we specify
‘variables’ as the vector to be integrated, and the other vectors as
<dependencies>. Here, for convenience, we introduce the variables sum1 and
sum2 which represent the two summations seen in equations 11. We also
introduce variables such as rhonm as a shorthand for rho(m => n, n => m).
Lastly we specify the variables we wish to write to memory in <output>. We are
only interested in the real part of $J_{mm}$ so we introduce a save variable
rhoReOut inside <moments>, and specify that our <dependencies> are the vector
<variables>.

###  1.1 Recreating the figures in the main text

In the GitHub repository [27] there are also folders for each of the figures,
allowing one to easily reproduce them. Figures 1, 2 and 3 consist of
Mathematica files with a single cell that needs to be evaluated. The Figure 4
directory consists of sub-directories for Figure4a and Figure 4b. Figure4a
contains the XMDS2 code and a Mathematica file for plotting the output. Simply
compile the program with xmds2 org-pol-lat.xmds and then run with ./org-pol-
lat. Figure4b contains an additional bash script which allows one to
automatically run the simulation several times, saving the output to a new
folder for each run. It can be run as ./stochastic.sh x where x is the number
of repeat simulations required. Figure5 contains very similar code to
Figure4a, except in this case M = 40 and we specify a single pump site with
pump = (m == M / 2) ? P0 : 0;

##  2 Checking the analytical calculations with QuantumAlegra.jl

To double check our analytical calculations we wrote a custom script with the
Julia package QuantumAlgebra.jl [28] (the package is maintained by Johannes
Feist from the Autonomous University of Madrid). The full script can be
downloaded [27], and could be helpful to anybody interested in calculating
mean-field equations of motions from Lindblad master equations.

###  2.1 Preliminaries

Consider that we have a known Hamiltonian that is a function of the bosonic
operators $a_{n}$. We have some arbitrary operators
$\hat{A}_{j}=f_{j}(a_{1},a_{2},\dots,a_{N})$, and we want to calculate the
mean-field dynamical equations for $A_{j}=\langle\hat{A}_{j}\rangle$. To
accomplish this we can use the identity
$\langle\hat{X}\rangle=\text{Tr}[\rho\hat{X}]$, which gives us

$\frac{\text{d}A_{j}}{\text{d}t}=\text{Tr}\left[\frac{\text{d}\rho}{\text{d}t}\hat{A}_{j}\right],$
(8)

as well as the Lindblad master equation

$\frac{\text{d}\rho}{\text{d}t}=\text{i}\left[\rho,H\right]+\sum_{n}\gamma_{n}\mathcal{L}(L_{n},\rho),$
(9)

where $L_{n}$ are the Lindblad jump operators of the system (in simple cases
we encounter just the operators $a_{n}$ and maybe $a_{n}^{\dagger}$), and

$\mathcal{L}(x,y)=\left(xyx^{\dagger}-\frac{1}{2}\left\\{x^{\dagger}x,y\right\\}\right).$
(10)

To calculate the dynamical equations we exploit the fact that the trace is
invariant under cyclic permutations. After substituting equation (2) into
equation (1), this allows us to rewrite the mean-field equations as

$\frac{\text{d}A_{j}}{\text{d}t}=\left\langle-\text{i}[\hat{A}_{j},H]+\sum_{n}\gamma_{n}\mathcal{M}(L_{n},\hat{A}_{j})\right\rangle$
(11)

where

$\mathcal{M}(x,y)=\left(x^{\dagger}yx-\frac{1}{2}\left\\{x^{\dagger}x,y\right\\}\right).$
(12)

###  2.2 The code

We transcribe the function $\mathcal{M}(x,y)$ to code as

⬇ function M(x,y)  adjoint(x)*y*x - 0.5*adjoint(x)*x*y - 0.5*y*adjoint(x)*x
end

and implement the RHS of equation 11 as

⬇ function ddt(H,A,$\gamma$,L)  -im*comm(A,H) + sum( map(($\gamma$,L) ->
$\gamma$ * M(L,A), $\gamma$, L))  end

where H is the Hamiltonian, A is the operator for which we want to calculate
the mean-field dynamics (e.g. the number of excitations adag(1)a(1)), and
$\gamma$ (resp. L) is a vector whose components are the $\gamma_{n}$ (resp.
$L_{n}$) we see in equation 11.

###  2.3 A simple example

Consider the Hamiltonian

$H=\omega_{1}a_{1}^{\dagger}a_{1}+\omega_{2}a_{2}^{\dagger}a_{2}+\Omega(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1}).$

The system experiences losses with rates $\gamma_{1}$ and $\gamma_{2}$, and
the 1st site is pumped incoherently with a strength $P$. Thus the master
equation can be written as

$\frac{\text{d}\rho}{\text{d}t}=\text{i}\left[\rho,H\right]+\gamma_{1}\mathcal{L}(a_{1},\rho)+\gamma_{2}\mathcal{L}(a_{2},\rho)+P\left[\mathcal{L}(a_{1},\rho)+\mathcal{L}(a_{1}^{\dagger},\rho)\right].$

The mean-field equation for the number of excitations on site 1 can be
calculated as

⬇ H = param(:$\omega_{1}$)*adag(1)*a(1) + param(:$\omega_{2}$)*adag(2)*a(2)  H
= H + param(:$\Omega$)*(adag(1)*a(2) + adag(2)*a(1))  $\gamma$ =
[param(:$\gamma_{1}$),param(:$\gamma_{2}$),param(:P)]  L =
[a(1),a(2),a(1)+adag(1)];  ddt(H,adag(1)a(1),$\gamma$,L)

##  3 Derivation from the exciton-photon basis

In the main text we consider a Hamiltonian that describes an excitonic
reservoir, polaritonic modes, and vibrons. This model could also describe
upper and lower polaritons, in the regime where the lower polaritons are
strongly photonic. In this case, the Hamiltonian in the main text can be
derived from the Hamiltonian in the exciton-photon basis
$H=H_{0}+H_{\mathrm{F}}+H_{\mathrm{I}}$, where

$\displaystyle H_{0}$
$\displaystyle=\sum_{m}\left(\omega_{\mathrm{cav}}c_{m}^{\dagger}c_{m}+\omega_{\mathrm{exc}}d_{m}^{\dagger}d_{m}+\omega_{\mathrm{V}}b_{m}^{\dagger}b_{m}\right)+\Omega_{\mathrm{R}}\sum_{m}\left(c_{m}^{\dagger}d_{m}+d_{m}^{\dagger}c_{m}\right),$
(13) $\displaystyle H_{\mathrm{F}}$
$\displaystyle=g_{0}\sum_{m}d_{m}^{\dagger}d_{m}\left(b_{m}+b_{m}^{\dagger}\right),$
(14) $\displaystyle H_{\mathrm{I}}$ $\displaystyle=\frac{1}{2}\sum_{m\neq
n}\Omega_{mn}^{0}\left(c_{m}^{\dagger}c_{n}+c_{n}^{\dagger}c_{m}\right).$ (15)

Here, $c^{\dagger}$ ($d^{\dagger}$) creates a photon (exciton) with frequency
$\omega_{\text{cav}}$ ($\omega_{\text{exc}}$), $\Omega_{\text{R}}$ is the
Rabi-frequency, $g_{0}$ is the interaction constant between excitons and
vibrons, and $\Omega_{mn}^{0}$ is the coupling between the $m$th and $n$th
cavity. The operators that diagonalise the single-cavity light-matter part of
the Hamiltonian are $\psi_{m}=d_{m}\cos\phi-c_{m}\sin\phi$ and
$a_{m}=c_{m}\cos\phi+d_{m}\sin\phi$, where
$\phi=(1/2)\mathrm{atan2}\left(\omega_{\mathrm{cav}}-\omega_{\mathrm{exc}},2\Omega_{\mathrm{R}}\right)$.
In this basis we have

$H_{0}=\sum_{m}\left(\omega_{\mathrm{C}}\psi_{m}^{\dagger}\psi_{m}+\omega_{\mathrm{R}}a_{m}^{\dagger}a_{m}+\omega_{\mathrm{V}}b_{m}^{\dagger}b_{m}\right),$
(16)

where

$\omega_{\mathrm{R},\mathrm{C}}=\frac{\omega_{\mathrm{exc}}+\omega_{\mathrm{cav}}}{2}\pm\sqrt{\left(\frac{\omega_{\mathrm{exc}}-\omega_{\mathrm{cav}}}{2}\right)^{2}+\Omega_{\mathrm{R}}^{2}}.$
(17)

The Fröhlich Hamiltonian becomes

$H_{\mathrm{F}}=g\sum_{m}\left(a_{m}b_{m}^{\dagger}\psi_{m}^{\dagger}+a_{m}^{\dagger}b_{m}\psi_{m}\right),$
(18)

where we have neglected nonresonant terms, and $g=g_{0}\sin\phi\cos\phi$. In
this work, we assume that $\cos^{2}\theta\ll\sin^{2}\theta$, in which case the
inter-microcavity interaction term becomes

$H_{\mathrm{I}}=\frac{1}{2}\sum_{m\neq
n}\Omega_{mn}\left(\psi_{m}^{\dagger}\psi_{n}+\psi_{n}^{\dagger}\psi_{m}\right),$
(19)

where $\Omega_{mn}=\Omega_{mn}^{0}\sin^{2}\phi$. With these approximations we
recover the Hamiltonian in the main text.