# Contrastive Multi-view Framework for Customer
Lifetime Value Prediction

Chuhan Wu1, Jingjie Li1, Qinglin Jia1, Hong Zhu2, Yuan Fang2, Ruiming Tang1
1Noah’s Ark Lab, Huawei 2Consumer Cloud Service Interactive Media BU, Huawei
wuchuhan1, lijingjie1, jiaqinglin2, zhuhong8, frank.fy<EMAIL_ADDRESS>
(2023)

###### Abstract.

Accurate customer lifetime value (LTV) prediction can help service providers
optimize their marketing policies in customer-centric applications. However,
the heavy sparsity of consumption events and the interference of data variance
and noise obstruct LTV estimation. Many existing LTV prediction methods
directly train a single-view LTV predictor on consumption samples, which may
yield inaccurate and even biased knowledge extraction. In this paper, we
propose a contrastive multi-view framework for LTV prediction, which is a
plug-and-play solution compatible with various backbone models. It synthesizes
multiple heterogeneous LTV regressors with complementary knowledge to improve
model robustness and captures sample relatedness via contrastive learning to
mitigate the dependency on data abundance. Concretely, we use a decomposed
scheme that converts the LTV prediction problem into a combination of
estimating consumption probability and payment amount. To alleviate the impact
of noisy data on model learning, we propose a multi-view framework that
jointly optimizes multiple types of regressors with diverse characteristics
and advantages to encode and fuse comprehensive knowledge. To fully exploit
the potential of limited training samples, we propose a hybrid contrastive
learning method to help capture the relatedness between samples in both
classification and regression tasks. We conduct extensive experiments on a
real-world game LTV prediction dataset and the results validate the
effectiveness of our method. We have deployed our solution online in Huawei’s
mobile game center and achieved 32.26% total payment amount gains.

Lifetime Value Prediction, Multi-view, Contrastive Learning

††copyright: rightsretained††journalyear: 2023††price: 15.00††doi:
10.1145/3534678.3539119††isbn: 978-1-4503-9385-0/22/08††ccs: Information
systems Personalization††ccs: Information systems Business intelligence††ccs:
Computing methodologies Multi-task learning

##  1\. Introduction

Customer lifetime value (LTV) reflects the expected revenue contributed by a
customer in a long-term relationship with a business (Jain and Singh, 2002;
Gupta et al., 2006). With the guidance of LTV in business planning and user
segmentation, service providers can make more informed marketing decisions and
provide better personalized services to increase user retention and reduce
churn (Berger and Nasr, 1998; Mani et al., 1999; Rosset et al., 2003;
Venkatesan and Kumar, 2004; Glady et al., 2009a; Pollak, 2021). Thus, accurate
LTV prediction is playing an increasingly important role in the digital
economy and intelligent management (Malthouse and Blattberg, 2005; Chamberlain
et al., 2017; Dahana et al., 2019; Win and Bo, 2020).

In many real-world applications, users’ consumption behaviors are sparse,
volatile, and noisy in their nature (Xing et al., 2021; Li et al., 2022). For
example, as shown in Fig. 1, in the App usage scenario, users may only choose
a few Apps from the various candidates recommended by the App store, and they
may order paid services in fewer Apps after trials. In our practice, only
about 2% of App downloads are finally converted into valid payments within a
month. Furthermore, users’ purchase behaviors may be affected by various
subject and objective factors, such as unplanned purchase intent and sales
promotion (Blattberg et al., 2009; Akram et al., 2018; Zhang et al., 2020). As
a result, we can observe a quite wide-range, long-tailed, and noisy LTV
distribution (see our dataset statistics for an intuitive example). The heavy
sparsity and noise of user consumption samples make LTV prediction a rather
challenging problem.

Figure 1. The consumption behavior patterns of a user in a typical App usage
scenario.

Considering the sparsity of samples with non-zero LTVs, many methods decompose
the LTV prediction problem into a binary classification task to predict
premium users and a regression task to estimate potential monetary values
(Wang et al., 2019). Some of them use two-stage strategies to model the two
tasks in a cascaded way and use machine learning models such as random forests
(Vanderveld et al., 2016) and XGBoost (Drachen et al., 2018) to process input
features. However, maintaining multiple models in these two-stage methods may
yield higher complexity and heavier error propagation problem. In recent
years, several deep learning-based methods provide end-to-end solutions to LTV
prediction. For example, ZILN (Wang et al., 2019) unifies the two tasks in a
single model via multi-task learning, where the final predicted LTV is given
by the multiplication between the purchase propensity and the expected
purchase amount. ODMN (Li et al., 2022) models the temporal dependencies
between LTVs, and uses a multi-distribution multi-expert module that divides
the imbalanced LTV distribution into several smoother sub-distributions and
assigns them to different experts. However, the single-view LTV regressors in
these methods may still be severely perturbed and even biased by the outlier
and noisy samples, especially when there are only very limited positive
samples for model optimization.

In this paper, we propose a Contrastive Multi-view framework for LTV
prediction named CMLTV, which can effectively mitigate the impact of the
consumption data sparsity, volatility, and noise on model learning. The core
of CMLTV is a multi-view regression framework, where multiple heterogeneous
regressors with diverse characteristics and advantages are jointly optimized
to extract complementary and robust knowledge. Their predictions are
synthesized into a unified score, which is further combined with the purchase
probability estimated by a purchase classifier to generate the final LTV. To
reduce the model’s dependency on data volume, we propose a hybrid contrastive
learning method to help capture the relatedness between samples in both
classification and regression tasks to fully exploit data potentials. By
jointly regularizing the purchase classifier and multi-view regressors, the
model can be aware of the inherent relatedness among samples and make more
accurate predictions. Our offline experiments are conducted on a real-world
game App user LTV prediction dataset collected from Huawei’s mobile game
center, which shows the effectiveness of CMLTV. Further online A/B test on
this platform also shows 33.26% total payment amount gains. Up to now, CMLTV
has been deeply involved to provide stable and high-quality game App
suggestion services for hundreds of millions of users.

The contributions of this paper are listed as follows:

  * •

We propose a multi-view LTV prediction framework that can encode complementary
knowledge via multiple heterogeneous regressors to confront data noise and
imbalance.

  * •

We propose a hybrid contrastive learning method that enables the model to
organically exploit the relatedness between samples in both classification and
regression tasks to mitigate the impact of data sparsity.

  * •

We conduct extensive experiments in both offline and online environments, and
the results fully validate the effectiveness of our method. Our solution has
been deployed online for personalized game App suggestions, serving hundreds
of millions of mobile users.

##  2\. Related Work

###  2.1. LTV Prediction

LTV prediction is a classic problem in business management (Kumar et al.,
2004). Pioneer research on customer LTV prediction focuses on building
probabilistic models based on observed data (Jasek et al., 2019). For example,
Pareto/NBD (Schmittlein et al., 1987) is a canonical model for customer
analysis based on historical transactions, which forecasts the future activity
and purchase frequency of users by modeling user behaviors with certain
stochastic processes. Based on this framework, Fader et al. (2005) studied
grouping users based on the recency, frequency, and monetary value of their
transactions, and then estimated the LTV for each customer cohort. These
probability-based methods require prior knowledge of user behaviors to build
basic model assumptions, which are highly sensitive to data noise and
volatility. In complex scenarios, it is also difficult to simulate user
behaviors via well-established stochastic processes.

Another major research line of LTV prediction is using machine learning
techniques to build data-driven models from historical logs (Tkachenko, 2015;
Chen et al., 2018; Sabbeh, 2018; del Río et al., 2019; Burelli, 2019; Jiang,
2020; Tekin et al., 2022). Considering the characteristics of LTV
distribution, several methods decompose the task into two stages, i.e., stage
one predicts the purchase propensity and stage two estimates the LTVs of users
who are predicted to purchase in stage one. For example, Vanderveld et al.
(2016) proposed to use random forests in both stages, which are built on
various features including user engagement measurements, user experience,
historical behaviors, and demographics. Drachen et al. (2018) proposed to use
XGBoost models to implement a game LTV prediction framework that involves
various gameplay, social, and purchase features. These two-stage methods need
to maintain multiple cascaded models, which may lead to high computational
costs and online latency.

In recent years, several deep learning-based methods unify the two stages into
a single model. For example, ZILN (Wang et al., 2019) uses a multi-task
learning framework that optimizes a binary classifier to predict purchase
propensity and a lognormal distribution-based regressor to predict the
expected payment amount. The final LTV is given by multiplying the predicted
purchase propensity and payment amount. ODMN (Li et al., 2022) further
considers the sequential dependencies between LTV in different periods. At
each prediction point, it uses a multi-distribution multi-expert module that
predicts the classification probabilities of different LTV ranges and uses
them to adaptively select proper experts to predict LTV in certain ranges.
Unfortunately, due to the noisy, imbalanced, and volatile nature of purchase
data, the single-view LTV regressors in these methods may still be heavily
affected and even biased, especially when there are insufficient samples to
draw a panorama of data distribution. Our work incorporates a multi-view
regression framework to extract complementary information and fuses it into
more robust and accurate knowledge. It is remarkable that several works in
recent years studied LTV prediction in other aspects, such as user
representation (Xing et al., 2021), purchase sequence modeling (Chen et al.,
2018; Bauer and Jannach, 2021), and social information mining (Zhao et al.,
2023). Since these scenario-specific methods bring some additional constraints
to data formats, they do not fall into the same research line as our work.

Figure 2. The framework of our CMLTV method.

###  2.2. Contrastive Learning in Recommendation

Contrastive learning is a widely employed technique in the recommendation
field (Wei et al., 2021; Yu et al., 2022), and here we briefly introduce a few
relevant works. A common contrastive paradigm in recommendation is comparing
each positive sample with one or more negative samples, e.g., BPR (Rendle et
al., 2009) and InfoNCE (Oord et al., 2018) losses. This paradigm is a popular
choice in learning-to-rank systems in various recommendation scenarios, such
as e-commerce (Xie et al., 2022), video (Yi et al., 2022), and news feed (Wu
et al., 2019). However, in online advertising applications, both the orders of
candidate ads and the absolute values of predicted factors (e.g., pCVR and
pLTV) matter in the bidding process (Guo et al., 2017; Pan et al., 2021).
Thus, the pair-wise contrastive learning paradigm may be inappropriate in
advertising scenarios due to its distortion effects on predicted scores. To
adopt the strong representation ability of contrastive learning in advertising
systems, Wang et al. (2022) proposed a CL4CTR framework that uses data
augmentation to construct contrast pairs. However, it cannot effectively
exploit the relatedness between different samples. Moreover, most studies of
contrastive learning in recommendation focus on classification problems (Lin
et al., 2022; Bai et al., 2022), which are not fully compatible with LTV
prediction that involves a regression task. We propose a hybrid contrastive
learning method that simultaneously guides model optimization in both
classification and regression task, which is especially suitable for LTV
prediction to alleviate the impact of data scarcity.

##  3\. Methodology

Here we introduce our CMLTV approach for customer LTV prediction. Since there
are several different formulations of LTV prediction in existing literature,
we first present the definition of the problem studied in our paper. We then
describe our model in detail and show how to optimize it to generate LTV
predictions.

###  3.1. Problem Definition

Denote a user and a certain product/service by $u$ and $i$, respectively.
Given the features of $u$ and $i$ (as well as some context features if
available), the goal of an LTV prediction model is to estimate the payment
amount of $u$ on $i$ within a certain time span $T$ (e.g., 30 days). If the
user $u$ does not purchase $i$ within $T$, the corresponding LTV in this
period is labeled as zero. If there are multiple transactions concerning $u$
and $i$ in this time span, the LTV label is the summation of their monetary
values. After model learning on historical transaction data, it is expected to
predict future LTVs as a reference for decision making and personalized
advertising.

###  3.2. CMLTV Framework

The overall framework of our CMLTV method is shown in Fig. 2. When a batch of
training samples arrives, we first use a multi-view LTV prediction model to
generate the purchase probability and multiple LTV regression results for each
sample, and then apply a hybrid contrastive learning strategy to samples in
this batch to capture their inherent relatedness. Their details are described
below.

####  3.2.1. Multi-View LTV Prediction

The multi-view LTV prediction module uses heterogeneous regressors111The MLP
networks in regressors may contain one or more layers, but we denote them as
“MLP” for consistency. with diverse characteristics to profile the input
sample in different aspects. Denote the input features of a sample by
$\mathbf{x}$.222We assume that discrete features have been converted into
embeddings. Firstly, a basic feature representation model is used to learn the
hidden representations $\mathbf{h}$ of input features by modeling their
interactions. Note that our framework does not restrict the architecture of
this backbone model, and it can be implemented by various off-the-shelf
structures such as MLP, DeepFM (Guo et al., 2017), DCN (Wang et al., 2017),
and DCNv2 (Wang et al., 2021).

Next, a purchase classifier is applied to the hidden representation
$\mathbf{h}$ to predict the purchase probability $\hat{p}$ as follows:

(1) $\displaystyle\mathbf{h}_{p}$
$\displaystyle=\rm{ReLU}(\mathbf{W}_{p}\mathbf{h}+\mathbf{b}_{p}),$
$\displaystyle\hat{p}$
$\displaystyle=\sigma(\mathbf{w}_{p}\mathbf{h}_{p}+b_{p}),$

where $\mathbf{W}_{p}$, $\mathbf{w}_{p}$, $\mathbf{b}_{p}$, and $b_{p}$ are
parameters, and $\sigma(\cdot)$ stands for the sigmoid function. The
probability output by the classifier is used to compute the cross-entropy loss
$\mathcal{L}_{p}$ as follows:

(2) $\mathcal{L}_{p}=-z\log\hat{p}-(1-z)\log{(1-\hat{p})},$

where $z$ is the binary purchase label of a sample.

We then estimate the payment amount of each sample in a multi-view manner.
Here we include three different types of regressors. The first one is a
distribution-based regressor, which aims to approximate real LTV distributions
with well-formulated probabilistic models. Although it can be implemented by
any probability distributions, we recommend using some canonical ones such as
lognormal and gamma distributions that have been shown effective in
characterizing long-tail property (Kim and Yum, 2008; Glady et al., 2009b;
Yang et al., 2022). Here we denote the normalized probability density function
of the selected distribution by $f(x)$. It is parameterized by a variable333It
can be either a scalar or a vector, depending on the parameters required by
$f(x)$. $\mathbf{\theta}$, which is learned from $\mathbf{h}$ as follows:

(3) $\displaystyle\mathbf{h}_{d}$
$\displaystyle={\rm{ReLU}}(\mathbf{W}_{d}\mathbf{h}+\mathbf{b}_{d}),$
$\displaystyle\mathbf{\theta}^{\prime}$
$\displaystyle=\mathbf{w}_{d}\mathbf{h}_{d}+b_{d},$
$\displaystyle\mathbf{\theta}$
$\displaystyle=\log[1+\exp(\mathbf{\theta}^{\prime})],$

where $\mathbf{W}_{d}$, $\mathbf{w}_{d}$, $\mathbf{b}_{d}$, and $b_{d}$ are
learnable parameters. The loss function $\mathcal{L}_{d}$ of this regressor is
the negative log-likelihood of the probability density given by $f(x)$, which
is formulated as follows:

(4) $\mathcal{L}_{d}=-\log f(y),$

where $y$ is the LTV label of a sample. By optimizing this loss function, the
model is encouraged to find better $\mathbf{\theta}$ to parameterize $f(x)$
that can maximize the LTV probability density of training samples. In the test
phase, it uses the expectation of the distribution $f(x)$ parameterized by the
$\mathbf{\theta}$ obtained from each sample as the prediction. In our
approach, we use gamma distribution to instantiate the function $f(x)$.
Therefore, the variable $\mathbf{\theta}$ is a vector including two elements,
i.e., a shape parameter and a rate parameter, and the predicted LTV score is
the division between them.

The second regressor is logarithm-based, which aims to predict the LTV score
in the log scale. Although customer LTVs can range from zero to million-level
numbers, they can be only low single digits in the logarithm view, which are
suitable for common neural networks to handle. Concretely, the output
$\hat{y}_{l}$ of this regressor is obtained from $\mathbf{h}$ through an MLP
module as follows:

(5) $\displaystyle\mathbf{h}_{l}$
$\displaystyle={\rm{ReLU}}(\mathbf{W}_{l}\mathbf{h}+\mathbf{b}_{l}),$
$\displaystyle\hat{y}_{l}$
$\displaystyle={\rm{ReLU}}(\mathbf{w}_{l}\mathbf{h}_{l}+b_{l}),$

where $\mathbf{W}_{l}$, $\mathbf{w}_{l}$, $\mathbf{b}_{l}$, and $b_{l}$ are
MLP parameters. The loss function $\mathcal{L}_{l}$ of this regressor is the
standard Mean Squared Error (MSE), which is computed as follows:

(6) $\mathcal{L}_{l}=(\hat{y}_{l}-\lg(1+y))^{2}.$

Here we add one to the label to ensure the non-negativity.

The third regressor is classification-based, which first converts the
regression task into a classification problem and then reconstructs a real-
valued score from the class probabilities as the final prediction. Since LTV
distribution is usually long-tailed, it divides the whole range into several
class bins using logarithm operations. We use base 2 logarithms for binning
rather than larger bases because it can keep more fine-grained numerical
information.444While it is possible to choose smaller log bases to achieve
finer-grained binning, it may lead to an excess of class numbers and even
redundant classes. More specifically, the corresponding class ID $c$ of LTV
label $y$ is obtained as follows:

(7) $c=\lfloor\log_{2}(1+y)\rfloor.$

The class probability vector $\mathbf{\hat{y}}_{c}$ is predicted as follows:

(8) $\displaystyle\mathbf{h}_{c}$
$\displaystyle={\rm{ReLU}}(\mathbf{W}_{c}\mathbf{h}+\mathbf{b}_{c}),$
$\displaystyle\mathbf{\hat{y}}_{c}$
$\displaystyle={\rm{Softmax}}(\mathbf{V}_{c}\mathbf{h}_{c}+\mathbf{v}_{c}),$

where $\mathbf{W}_{c}$, $\mathbf{V}_{c}$, $\mathbf{b}_{c}$, and
$\mathbf{v}_{c}$ are classifier parameters. The loss function
$\mathcal{L}_{c}$ for model learning is the following multi-class cross-
entropy:

(9) $\mathcal{L}_{c}=-\sum_{i=1}^{C}c_{i}\log\mathbf{\hat{y}}_{c}^{i},$

where $C$ is the class number, $c_{i}$ and $\hat{y}_{c}^{i}$ is the real and
predicted label for the $i$-th class, respectively. The final real-valued LTV
prediction $\hat{y}^{\prime}_{c}$ is reconstructed from the probability vector
$\mathbf{\hat{y}}_{c}$ as follows:

(10)
$\hat{y}^{\prime}_{c}=\sum_{i=1}^{C}\frac{(2^{i}-1+2^{i+1}-2)}{2}\mathbf{\hat{y}}_{c}^{i}=\sum_{i=1}^{C}\frac{(2^{i}+2^{i+1}-3)}{2}\mathbf{\hat{y}}_{c}^{i}.$

In this way, we use the expectation as LTV prediction, where the center of
each bin is used to represent its expected LTV.

####  3.2.2. Hybrid Contrastive Learning

Next, we introduce the hybrid contrastive learning mechanism in our approach,
which aims to build organic connections between samples in the same batch to
better exploit the potential of limited training data. Assume there are $K$
samples within a batch. Their predicted purchase probabilities and three types
of LTV scores are denoted by $[\hat{p}_{1},...,\hat{p}_{K}]$,
$[\hat{y}_{d,1},...,\hat{y}_{d,K}]$, $[\hat{y}_{l,1},...,\hat{y}_{l,K}]$,
$[\hat{y}_{c,1},...,\hat{y}_{c,K}]$, respectively. Inspired by BPR loss, we
encourage the purchase probabilities of positive samples to be larger than
negative samples. However, it is infeasible to directly compare each pair of
samples because the noisy samples may mislead the model to generate low-
quality contrastive pairs, which are harmful to model optimization.
Fortunately, although parts of the samples are uninformative, positive samples
should be assigned higher purchase propensities than negative samples on
average. Thus, we first compute the average purchase probabilities of positive
and negative samples, which are denoted by $\hat{p}_{+}$ and $\hat{p}_{-}$,
respectively. We design a classification contrastive loss $\mathcal{L}^{c}$ by
comparing $\hat{p}_{+}$ and $\hat{p}_{-}$ as follows:

(11)
$\mathcal{L}^{c}=-\log\sigma[\sigma^{-1}(\hat{p}_{+})-\sigma^{-1}(\hat{p}_{-})],$

where $\sigma^{-1}$ means the inverse function of sigmoid.

Motivated by the findings in prior work (Borle et al., 2008), we assume that
users with higher LTVs may also have higher purchase probabilities. This
assumption is probably true in various scenarios because large consumption
values are usually caused by multiple purchase behaviors. Thus, we propose a
regression contrastive learning method to regularize the regression results to
be positively correlated with the predicted purchase probabilities. We take
the scores output by the distribution-based regressor as an example to
elaborate this process. Its corresponding regression contrastive loss
$\mathcal{L}^{r}_{d}$ is computed as follows:

(12)
$\mathcal{L}^{r}_{d}=-\frac{1}{K^{2}}\sum_{i=1}^{K}\sum_{j=1}^{K}(\hat{p}_{i}-\hat{p}_{j})[\lg(1+\hat{y}_{d,i})-\lg(1+\hat{y}_{d,j})].$

Here we use base 10 logarithms to control the scale of contrastive loss.
Similarly, we compute the regression contrastive losses $\mathcal{L}^{r}_{l}$
and $\mathcal{L}^{r}_{c}$ based on the other two types of regression results.
In this way, the relatedness between different samples in the same batch is
explicitly encoded into model learning, meanwhile the classification and
regression parts are naturally connected by these losses so that their encoded
knowledge can be exchanged and shared.

###  3.3. Model Training and Prediction

Finally, we introduce the details of model training and testing in our
framework. Our method uses a multi-task learning framework to unify different
types of loss functions. The overall loss $\mathcal{L}$ for model training is
a combination of the binary classification loss, three types of regression
losses, contrastive classification loss, and three types of contrastive
regression losses. Note that the binary classification loss and the
contrastive losses are computed on all samples, while the regression losses
are only activated on positive samples. This is because the heavy data
imbalance shall damage the accuracy of regression models. Therefore, the loss
$\mathcal{L}$ is formulated as follows:

(13)
$\mathcal{L}=\mathcal{L}^{c}+\mathcal{L}^{r}_{d}+\mathcal{L}^{r}_{l}+\mathcal{L}^{r}_{c}+\sum_{i}[\mathcal{L}_{p}(i)]+\sum_{y_{i}>0}[\mathcal{L}_{d}(i)+\mathcal{L}_{l}(i)+\mathcal{L}_{c}(i)],$

where $\mathcal{L}_{p}(i)$, $\mathcal{L}_{d}(i)$, $\mathcal{L}_{l}(i)$, and
$\mathcal{L}_{c}(i)$ represent the corresponding loss on the $i$-th sample,
and $y_{i}$ is its LTV label. By optimizing this objective function, the model
can fuse the knowledge extracted by the multi-view framework to generate
accurate predictions.

After model convergence, we use the multi-view LTV prediction part for
inference. Motivated by the ZILN (Wang et al., 2019) framework, we multiply
the purchase probability and the regression scores. The final LTV score
$\hat{y}$ is generated as follows:

(14)
$\hat{y}=\hat{p}\cdot[\alpha\hat{y}_{d}+\beta\hat{y}_{l}+(1-\alpha-\beta)\hat{y}_{c}],$

where $\alpha$ and $\beta$ are weighting coefficients.555We do not fuse the
three scores in the model training phase to learn these coefficients because
it will lead to performance degradation.

###  3.4. Model Complexity Analysis

Finally, we provide some analysis of the theoretical computational complexity
of our approach. The complexity of the multi-view LTV prediction part mainly
depends on the architecture of the basic model. If it is implemented by an MLP
model, the computational cost of this module is $O(Kd^{2})$, where $d$ is the
hidden dimension. The computational cost of our hybrid contrastive learning
part is mainly brought by computing the regression contrastive loss, whose
complexity is $O(K^{2})$. Thus, the overall computational complexity of our
framework is $O(Kd^{2}+K^{2})$. If the batch size $K$ is much larger than the
hidden dimension $d$, then the bottleneck is the contrastive learning
part.666In our experiments, the contrastive learning mechanism actually
occupies around 58% of the training time. Therefore, it is preferable to set a
moderate batch size due to efficiency considerations.

##  4\. Experiments

Figure 3. The distribution of non-zero LTVs in our dataset.

###  4.1. Dataset

Although there are a few public datasets for LTV prediction, they may not well
reflect the characteristics of many scenarios where positive samples are much
sparser than negatives.777Some e-commerce datasets only contain users’
transaction histories (i.e., positive samples), where negative samples are not
included. Thus, we collect a dataset from Huawei’s mobile game center for
offline experiments. It is sampled from the aggregated App consumption logs
from Aug. 1st to Aug. 31st, 2022, and it contains 642k positive samples and
30.1m negative samples in total. The LTV labels are counted within the 30 days
after download. For example, the LTV label of a sample generated on Aug. 1st
is the accumulated consumption amount from Jul. 2nd to Aug. 1st. Each sample
is associated with several numerical features (e.g. historical LTVs),
categorical features (e.g., App IDs and categories), and binary features
(e.g., user segments). The distribution of non-zero LTVs (the values are
normalized in the log scale) is illustrated in Fig. 3. We find it
approximately obeys a long-tail distribution but there are multiple peaks on
the curve, which is possibly caused by the institutionalized price levels on
mobile App stores. The samples on the last day are used for test and we
randomly sample 10% of training data for validation.

Table 1. LTV prediction performance of different methods on all samples or
positive samples only.

Methods | All Samples |  Positive Samples
---|---|---
| RMSE | MAE | Pearson | Spearman | R2_score | AUC | RMSE | MAE | Pearson | Spearman | R2_score
Linear | 247.16 | 9.809 | 0.1422 | 0.1194 | 0.01958 | 0.7517 | 1770.4 | 236.85 | 0.3229 | 0.2927 | 0.02588
MLP | 248.11 | 20.42 | 0.1257 | 0.0740 | 0.01209 | 0.6549 | 1789.0 | 232.96 | 0.3249 | 0.4772 | 0.00538
RF | 248.87 | 4.774 | 0.0949 | 0.1352 | 0.00599 | 0.7314 | 1803.0 | 241.59 | 0.2492 | 0.0949 | 0.01035
XGBoost | 246.40 | 4.681 | 0.2548 | 0.1581 | 0.02565 | 0.8336 | 1786.3 | 241.43 | 0.2664 | 0.2424 | 0.00835
ZILN | 242.56 | 8.090 | 0.2421 | 0.1708 | 0.05576 | 0.8601 | 1713.0 | 241.37 | 0.3456 | 0.4345 | 0.07721
MDME | 243.71 | 7.969 | 0.2346 | 0.1668 | 0.04395 | 0.8482 | 1743.6 | 243.59 | 0.3286 | 0.4299 | 0.07843
CMLTV | 240.16 | 6.889 | 0.2495 | 0.1724 | 0.07436 | 0.8629 | 1710.5 | 235.74 | 0.3746 | 0.4382 | 0.10291
|  |  |  |  |  |  |  |  |  |  |

Table 2. Comparison of different frameworks in empowering different models.

Methods | All Samples |  Positive Samples
---|---|---
| RMSE | MAE | Pearson | Spearman | R2_score | AUC | RMSE | MAE | Pearson | Spearman | R2_score
MLP+ZILN | 242.56 | 8.090 | 0.2421 | 0.1708 | 0.05576 | 0.8601 | 1713.0 | 241.37 | 0.3456 | 0.4345 | 0.07721
MLP+MDME | 243.71 | 7.969 | 0.2346 | 0.1668 | 0.04395 | 0.8482 | 1743.6 | 243.59 | 0.3286 | 0.4299 | 0.07843
MLP+CMLTV | 240.16 | 6.889 | 0.2495 | 0.1724 | 0.07436 | 0.8629 | 1710.5 | 235.74 | 0.3746 | 0.4382 | 0.10291
DCN+ZILN | 242.14 | 7.932 | 0.2425 | 0.1713 | 0.05685 | 0.8614 | 1704.2 | 239.89 | 0.3467 | 0.4348 | 0.07895
DCN+MDME | 243.53 | 7.885 | 0.2367 | 0.1689 | 0.04647 | 0.8525 | 1739.8 | 242.13 | 0.3305 | 0.4304 | 0.07822
DCN+CMLTV | 239.95 | 6.842 | 0.2510 | 0.1740 | 0.07569 | 0.8643 | 1699.6 | 234.36 | 0.3757 | 0.4390 | 0.10441
AutoInt+ZILN | 242.34 | 7.994 | 0.2420 | 0.1711 | 0.05664 | 0.8615 | 1706.7 | 241.01 | 0.3449 | 0.4340 | 0.07599
AutoInt+MDME | 243.69 | 7.910 | 0.2342 | 0.1672 | 0.04592 | 0.8509 | 1741.4 | 242.96 | 0.3278 | 0.4284 | 0.07644
AutoInt+CMLTV | 240.02 | 6.856 | 0.2494 | 0.1730 | 0.07498 | 0.8633 | 1702.3 | 234.98 | 0.3751 | 0.4377 | 0.10284
|  |  |  |  |  |  |  |  |  |  |

###  4.2. Experimental Settings

In our experiments, we use a two-layer MLP with ReLU activation functions as
the basic feature representation model in all deep learning methods, if not
specified. Each layer in the backbone model is followed by a batch
normalization (Ioffe and Szegedy, 2015) operation. The hidden units of this
model are 256 and 128, respectively. The hidden dimension in the classifier
and regressors is 64. We use Adam (Bengio and LeCun, 2015) as the model
optimizer and its learning rate is 1e-3. The batch size is set to 10,000. We
use the early stopping mechanism in model training and its patient is 3
epochs. The coefficients $\alpha$ and $\beta$ are set to 0.3. These
hyperparameters are tuned on the validation set.

Although many previous works (Wang et al., 2019; Li et al., 2022) use Rooted
Mean Squared Error (RMSE) and Mean Absolute Error (MAE)888We do not use
normalized RMSE and MAE because we find the data is so imbalanced that even
normalized RMSE and MAE are larger than 1. Thus, we directly use the original
non-normalized values.  as the performance metric, in our practice we find
they are insufficient to comprehensively reflect prediction quality,
especially when negative samples are dominant. Thus, in our experiments we use
six metrics to comprehensively evaluate model performance, including RMSE,
MAE, Pearson correlation coefficient (denoted by Pearson), Spearman
correlation coefficient (denoted by Spearman), coefficient of determination
(denoted by R2_score), and Area Under ROC Curve (AUC).999The AUC score is
evaluated by the binary purchase label. To fully understand model performance
on different types of samples, we compute all metrics except AUC on all
samples or positive samples only, respectively. This is because the errors on
positive samples and negative samples may yield different volumes of loss in
practical scenarios. We repeat each experiment five times to mitigate
randomness and report the average performance.

###  4.3. Performance Evaluation

In our offline evaluation, we compare our CMLTV approach with the following
baselines:

  * •

Linear, which uses a multivariate linear regression model for LTV prediction;

  * •

MLP, which uses a multi-layer perceptron network for LTV prediction;

  * •

RF (Vanderveld et al., 2016), which uses two independent random forest models
for classification and regression;

  * •

XGBoost (Drachen et al., 2018), replacing the random forest with the XGboost
model in the above two-stage method;

  * •

ZILN (Wang et al., 2019), a multi-task learning method that unified binary
classification and lognormal distribution-based regression;

  * •

MDME (Li et al., 2022), a multi-distribution multi-expert method for LTV
prediction, which is the core module in ODMN (Li et al., 2022).

The performance of compared methods is listed in Table 1, from which we have
the following observations:

Figure 4. The performance of CMLTV and its variants without different
regression views. Figure 5. The performance of using different types of
regressors alone.

First, we find different methods may have different ranks in terms of
different metrics. For example, on all samples RF achieves a low MAE score but
performs poorly in terms of AUC. The mismatch between different metrics is
mainly due to the high data imbalance and the perturbation of outlier data.
Thus, it is insufficient to merely use RMSE and MAE to measure model
performance. Only if a model performs better than other models in most
metrics, we can be confident about its effectiveness. According to the long-
term experience in our practical scenarios, R2_score usually has the best
correlation with real online performance. RMSE and MAE are less informative
since they are highly sensitive to outliers. Industrial practitioners can
choose it as the primary metric if using multiple metrics is not supported.
Second, the Linear and MLP models perform poorly in many metrics. This is
because they are heavily perturbed and even biased by the extremely imbalanced
data. It seems that MLP is more easily to be affected by data imbalance than
the simple linear regression model, which may be because of the high non-
linearity of deep networks. This shows the rationality of decomposing the LTV
prediction problem into a classification task and a regression task in deep
learning-based methods. Third, single-stage deep learning methods such as ZILN
and MDME usually outperform two-stage methods (i.e., RF and XGBoost). This is
because the separate models in the two stages cannot share their mutual
knowledge and the errors produced by the first model are completely inherited
by the second one. Thus, unifying the two processes in a single model is a
better choice. Forth, our CMLTV method outperforms other baselines (at the
significant level of $p<0.05$ in two-sided t-test, evaluated on the best-
performed metrics). This may be because our approach better defends data noise
and volatility meanwhile being aware of samples’ relationships to overcome
data sparsity. Thus, our method is more robust and accurate than baselines.

Since ZILN, MDME, and CMLTV are independent on the base model, we use them to
empower several widely used base models, including MLP, DCN (Wang et al.,
2017), and AutoInt (song2019autoint). The results are shown in Table 2. We
observe that CMLTV consistently outperforms ZILN and MDME, and the DCN based
model achieves the best results. This further verifies the effectiveness and
generality of our approach in leveraging existing models that are well crafted
for other tasks such as CTR prediction in LTV prediction.

###  4.4. Effectiveness of Regression Views

In this section we verify the effectiveness of different regression views in
our multi-view LTV prediction module. We compare the results of CMLTV with its
ablations with different regression views removed, as shown in Fig. 4. In our
later experiments, we only report the results on all samples due to space
limitations. We find it interesting that different regressors have different
advantages in view of different metrics. For example, the AUC substantially
declines when the gamma distribution-based regressor is removed, while the
Pearson correlation coefficient degrades the most when the classification-
based regressor is removed. This phenomenon confirms that different regression
views indeed encode heterogeneous and complementary knowledge, and none of
them is redundant. Therefore, combining the predictions from multi-view
regressors can keep satisfactory performance in all metrics.

We then report the results of using different types of regressors alone in our
framework to support the choice of regressors in our approach. Here we compare
five different regressors, including: (1) lognormal distribution-based
regressor in ZILN; (2) gamma distribution-based regressor in our method; (3)
negative exponential distribution-based regressor; (4) logarithm-based
regressor in our method; (5) classification-based regressor in our method, and
the results are shown in Fig. 5. We find gamma distribution is more suitable
than lognormal and exponential distributions for LTV modeling in our scenario,
which is probably because the shape of distribution curve shown in Fig. 3
better matches the shape of gamma distribution. Thus, we select gamma
distribution-based regressor in our method. In addition, logarithm-based and
classification-based regressors show advantages and disadvantages in different
metrics, thereby they can be complementary when being optimized jointly. Based
on the above analysis, we combine these regressors in our multi-view framework
to comprehensively empower our model.

###  4.5. Effectiveness of Hybrid Contrastive Learning

In this section we analyze the influence of the two types of contrastive
learning losses in our approach. The results of using different contrastive
loss combinations in our method are shown in Fig. 6. We observe that both
losses are beneficial for model learning, while they show different effects on
model performance. The classification contrastive learning loss mainly
improves AUC, which is intuitive because it is applied to the binary
classification task. On the contrary, the regression contrastive loss can
substantially improve both AUC and other regression metrics. This is because
it builds connections between predicted purchase probabilities and LTV scores,
thereby both classification and regression tasks can benefit from the mutual
knowledge they encode. Moreover, combining both tasks in our method can
consistently achieve better results, which validates the effectiveness of our
hybrid contrastive learning method.

Figure 6. Effectiveness of the two contrastive learning losses in CMLTV.
Figure 7. Comparisons among different types of contrastive learning paradigms.

(a) Purchase probabilities (y-axis is in log-scale). (b) LTV scores (x-axis is
in log-scale).

Figure 8. The distribution of predicted purchase probabilities and LTV scores.

We then compare our proposed contrastive learning approach with other mature
contrastive learning strategies, including the standard pair-wise method and
the data augmentation-based method introduced by (Wang et al., 2022). From the
results shown in Fig. 7, we find that the direct pair-wise contrast mechanism
has a negative effect on model performance. This is intuitive because it may
suffer from heavy data noise and unwanted prediction bias brought by the
contrast. The augmentation-based contrast method is also inferior to our
hybrid contrast method. This is because the former can only enhance the
discriminativeness of hidden representations and is agnostic to the
relatedness between real-valued labels of different samples. Our hybrid method
can collaboratively exploit the supervision signals across different samples
in both classification and regression tasks, which is especially powerful in
LTV prediction.

(a) All samples, linear scale. (b) Positive samples, linear scale. (c) All
samples, log scale. (d) Positive samples, log scale.

Figure 9. The Pearson correlation coefficients between different types of
scores in linear/log scales computed on all samples or positive samples only.
Figure 10. Influence of the prediction fusing coefficients on final
performance.

###  4.6. Prediction Results Analysis

We then present some analysis of scores predicted by CMLTV. We first
illustrate the distributions of purchase probabilities and the predicted LTV
scores in Fig. 8. Note that the LTV predictions that are smaller than 1 are
omitted. We find the predicted purchase probabilities of most samples are low.
This is intuitive because purchase behaviors are usually sparse. In addition,
we find the regression results given by different regressors have some
differences in their distributions, while their shapes are generally long-
tailed. Thus, by fusing the predictions from multiple heterogeneous
regressors, the model can have a more comprehensive and less biased
understanding of data distribution.

We also conduct further analysis of the prediction scores output by our
method. We compute the Pearson correlation coefficients between different
types of scores on all samples or positive samples only, where the LTV
regression scores are presented in both linear and log scales, as shown in
Fig. 9. We find that the regression results generated by different regressors
are highly correlated (especially on positive samples) but not identical. This
result shows that different regression modules have high agreements on most
samples but still keep their characteristics. Therefore, it is safe to
synthesize the predictions given by the three regressors. In addition, the
predicted purchase probabilities show strong correlations with predicted LTV
scores, especially in the log scale. This is mainly due to the regularization
of our contrastive regression loss so that predicted purchase probabilities
are positively correlated with LTV regressions.

###  4.7. Hyperparameter Analysis

In this section we study the influence of two key hyperparameters on our
approach, i.e., the prediction fusing coefficients $\alpha$ and $\beta$, which
controls the relative importance of distribution-based and logarithm-based
regressors in the final predictions. We vary the values of both coefficients
and the corresponding results are shown in Fig. 10. From the results, we find
the optimal points on the heatmaps of different metrics are diverse. For
example, the best MAE is achieved when $\alpha=1$, while the best Spearman
score is achieved when $\alpha=\beta=0$. Nonetheless, to achieve a good
tradeoff between different metrics, we prefer to select moderate values for
these coefficients. In our experiments we set $\alpha=\beta=0.3$, which
ensures that all three types of regressors can effectively contribute to the
final prediction meanwhile all types of metrics are satisfactory at this
point. In real-world scenarios, we recommend practitioners search these
hyperparameters according to the key metrics in the corresponding
applications.

Figure 11. The application scenario of our solution.

##  5\. Online Evaluation

We deployed our CMLTV method in Huawei’s mobile game center for game player
LTV prediction. The predictions given by our model are used to generate
personalized game suggestions on the “editor’s recommendations” game display
board, which is shown in Fig. 11. We performed a 70-day online A/B test on our
platform to verify the effectiveness of our algorithm.

###  5.1. Online Experimental Settings

The online A/B test lasts 70 days from Nov. 11, 2022, to Jan. 20, 2023. The
competitor of our model is a well-crafted baseline model, which is also
learned in an end-to-end manner via a multi-task learning framework. For fair
comparison, the baseline model and our approach share the same input features
and basic model backbones. Both models are retrained periodically on a dataset
aggregated from the transactions in the past month. Our model is deployed in a
virtual Linux computing node with an NVIDIA Tesla V100 GPU (32GB of memory).
In the online A/B test, 20% of randomly selected users are used as the
experimental group while another 20% of randomly selected users are reserved
as the control group. The traffics from the experimental and control groups
are served by our method and the baseline model, respectively. The predicted
LTVs are multiplied with the estimated CVR (in terms of App download) as the
ranking criteria to generate the final recommendation lists.

###  5.2. Online Results

Here we introduce the results of online A/B test. The traffics served by our
model achieved 32.26% payment amount improvements over the control group
traffics with a cost of 9% additional inference latency, which is fully
tolerable in our system.101010The scale of payment amount increase is
relatively larger than other types of common metrics such as CTR and CVR. This
is mainly because high consumption users could significantly boost the total
payment amount. In our further analysis, this significant improvement is
mainly due to the increased payment amount of users at low and medium
consumption levels (the revenue contributed by high consumption users is
similar for baseline and our approach). This finding indicates a critical
advantage of our approach in boosting the activeness of a broad range of
users, rather than only targeting high-spending users.

Due to the excellent performance of our CMLTV solution, it has become a core
module in the user-centric ecosystem of our App store, serving the main
traffic of hundreds of millions of users with various types of mobile devices.
Since both the multi-view regression and hybrid contrastive learning modules
in our approach are plug-and-play techniques, our solution has the potential
to empower various LTV prediction models in different scenarios. Compared with
many existing solutions, our approach has minimal requirements for the data
format and model pipeline, thereby can be adopted by different industrial
systems with light efforts, and we are currently working on transferring our
solution to other products such as online education and advertising.

##  6\. Conclusions and Limitations

In this paper, we present an industrial solution to customer lifetime value
prediction, named CMLTV. Different from conventional LTV prediction paradigms,
we propose a multi-view regression framework that assembles multiple
heterogeneous regression networks to capture complementary knowledge, which
endows the model with better prediction accuracy and robustness when handling
noisy and volatile consumption data. In addition, we propose a hybrid
contrastive learning mechanism to explicitly encode the relation information
between samples into model learning meanwhile building organic connections
between the binary purchase classifier and LTV regressors in our model. Hence,
the potential of training samples can be better exploited to help combat data
sparsity. Extensive experiments in both offline and online environments
validate the superiority of our approach over baseline methods. We have
deployed our solution online to serve hundreds of millions of our mobile users
and improve their experience.

Despite the effectiveness of our method, it still has the following
limitations. First, the model may underestimate the real LTV of some users
since the predicted purchase probabilities are usually always smaller than
one. Second, it is difficult for our method to predict very high LTVs, since
these outlier samples are very sparse in the training data. Third, our hybrid
contrastive learning method requires a large batch size, otherwise single
batch has too few positive samples and the loss is volatile. Finally, the LTVs
for cold-start users and items are hard to predict due to the lack of prior
knowledge (Singh et al., 2009). We will explore these directions in the
future.

## References

  * (1)
  * Akram et al. (2018) Umair Akram, Peng Hui, Muhammad Kaleem Khan, Yasir Tanveer, Khalid Mehmood, and Wasim Ahmad. 2018\.  How website quality affects online impulse buying: Moderating effects of sales promotion and credit card use.  _Asia Pacific Journal of Marketing and Logistics_ (2018).
  * Bai et al. (2022) Ting Bai, Yudong Xiao, Bin Wu, Guojun Yang, Hongyong Yu, and Jian-Yun Nie. 2022\.  A Contrastive Sharing Model for Multi-Task Recommendation. In _WWW_. 3239–3247.
  * Bauer and Jannach (2021) Josef Bauer and Dietmar Jannach. 2021.  Improved Customer Lifetime Value Prediction With Sequence-To-Sequence Learning and Feature-Based Models.  _TKDD_ 15, 5 (2021), 1–37.
  * Bengio and LeCun (2015) Yoshua Bengio and Yann LeCun. 2015.  Adam: A Method for Stochastic Optimization. In _ICLR_.
  * Berger and Nasr (1998) Paul D Berger and Nada I Nasr. 1998.  Customer lifetime value: Marketing models and applications.  _Journal of interactive marketing_ 12, 1 (1998), 17–30.
  * Blattberg et al. (2009) Robert C Blattberg, Edward C Malthouse, and Scott A Neslin. 2009.  Customer lifetime value: Empirical generalizations and some conceptual questions.  _Journal of Interactive Marketing_ 23, 2 (2009), 157–168.
  * Borle et al. (2008) Sharad Borle, Siddharth S Singh, and Dipak C Jain. 2008\.  Customer lifetime value measurement.  _Management science_ 54, 1 (2008), 100–112.
  * Burelli (2019) Paolo Burelli. 2019\.  Predicting customer lifetime value in free-to-play games.  In _Data analytics applications in gaming and entertainment_. Auerbach Publications, 79–107.
  * Chamberlain et al. (2017) Benjamin Paul Chamberlain, Angelo Cardoso, CH Bryan Liu, Roberto Pagliari, and Marc Peter Deisenroth. 2017.  Customer lifetime value prediction using embeddings. In _KDD_. 1753–1762.
  * Chen et al. (2018) Pei Pei Chen, Anna Guitart, Ana Fernández del Río, and Africa Periánez. 2018.  Customer lifetime value in video games using deep learning and parametric models. In _Big Data_. IEEE, 2134–2140.
  * Dahana et al. (2019) Wirawan Dony Dahana, Yukihiro Miwa, and Makoto Morisada. 2019\.  Linking lifestyle to customer lifetime value: An exploratory study in an online fashion retail market.  _Journal of Business Research_ 99 (2019), 319–331.
  * del Río et al. (2019) Ana Fernández del Río, Pei Pei Chen, and Africa Periánez. 2019\.  Profiling players with engagement predictions. In _CoG_. IEEE, 1–4.
  * Drachen et al. (2018) Anders Drachen, Mari Pastor, Aron Liu, Dylan Jack Fontaine, Yuan Chang, Julian Runge, Rafet Sifa, and Diego Klabjan. 2018\.  To be or not to be…social: incorporating simple social features in mobile game customer lifetime value predictions. In _ACSW_. ACM, 40:1–40:10.
  * Fader et al. (2005) Peter S Fader, Bruce GS Hardie, and Ka Lok Lee. 2005\.  RFM and CLV: Using iso-value curves for customer base analysis.  _Journal of marketing research_ 42, 4 (2005), 415–430.
  * Glady et al. (2009a) Nicolas Glady, Bart Baesens, and Christophe Croux. 2009a.  Modeling churn using customer lifetime value.  _European Journal of Operational Research_ 197, 1 (2009), 402–411.
  * Glady et al. (2009b) Nicolas Glady, Bart Baesens, and Christophe Croux. 2009b.  A modified Pareto/NBD approach for predicting customer lifetime value.  _Expert Systems with Applications_ 36, 2 (2009), 2062–2071.
  * Guo et al. (2017) Huifeng Guo, Ruiming Tang, Yunming Ye, Zhenguo Li, and Xiuqiang He. 2017.  DeepFM: a factorization-machine based neural network for CTR prediction. In _IJCAI_. 1725–1731.
  * Gupta et al. (2006) Sunil Gupta, Dominique Hanssens, Bruce Hardie, Wiliam Kahn, V Kumar, Nathaniel Lin, Nalini Ravishanker, and S Sriram. 2006\.  Modeling customer lifetime value.  _Journal of service research_ 9, 2 (2006), 139–155.
  * Ioffe and Szegedy (2015) Sergey Ioffe and Christian Szegedy. 2015.  Batch normalization: Accelerating deep network training by reducing internal covariate shift. In _ICML_. pmlr, 448–456.
  * Jain and Singh (2002) Dipak Jain and Siddhartha S Singh. 2002.  Customer lifetime value research in marketing: A review and future directions.  _Journal of interactive marketing_ 16, 2 (2002), 34–46.
  * Jasek et al. (2019) Pavel Jasek, Lenka Vrana, Lucie Sperkova, Zdenek Smutny, and Marek Kobulsky. 2019.  Comparative analysis of selected probabilistic customer lifetime value models in online shopping.  _Journal of Business Economics and Management_ 20, 3 (2019), 398–423.
  * Jiang (2020) Junxiang Jiang. 2020\.  A Study of Game Payment Data Mining: Predicting High-Value Users for MMORPGs. In _PAKDD Workshops_. Springer, 181–192.
  * Kim and Yum (2008) Jin Seon Kim and Bong-Jin Yum. 2008.  Selection between Weibull and lognormal distributions: A comparative simulation study.  _Computational Statistics & Data Analysis_ 53, 2 (2008), 477–485.
  * Kumar et al. (2004) V Kumar, Girish Ramani, and Timothy Bohling. 2004.  Customer lifetime value approaches and best practice applications.  _Journal of interactive Marketing_ 18, 3 (2004), 60–72.
  * Li et al. (2022) Kunpeng Li, Guangcui Shao, Naijun Yang, Xiao Fang, and Yang Song. 2022.  Billion-user Customer Lifetime Value Prediction: An Industrial-scale Solution from Kuaishou. In _CIKM_. ACM, 3243–3251.
  * Lin et al. (2022) Zihan Lin, Xuanhua Yang, Shaoguo Liu, Xiaoyu Peng, Wayne Xin Zhao, Liang Wang, and Bo Zheng. 2022.  Personalized Inter-Task Contrastive Learning for CTR&CVR Joint Estimation.  _arXiv preprint arXiv:2208.13442_ (2022).
  * Malthouse and Blattberg (2005) Edward C Malthouse and Robert C Blattberg. 2005.  Can we predict customer lifetime value?  _Journal of interactive marketing_ 19, 1 (2005), 2–16.
  * Mani et al. (1999) DR Mani, James Drew, Andrew Betz, and Piew Datta. 1999\.  Statistics and data mining techniques for lifetime value modeling. In _KDD_. 94–103.
  * Oord et al. (2018) Aaron van den Oord, Yazhe Li, and Oriol Vinyals. 2018.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ (2018).
  * Pan et al. (2021) Yujie Pan, Jiangchao Yao, Bo Han, Kunyang Jia, Ya Zhang, and Hongxia Yang. 2021\.  Click-through Rate Prediction with Auto-Quantized Contrastive Learning.  _arXiv preprint arXiv:2109.13921_ (2021).
  * Pollak (2021) Ziv Pollak. 2021\.  Predicting Customer Lifetime Values–ecommerce use case.  _arXiv preprint arXiv:2102.05771_ (2021).
  * Rendle et al. (2009) Steffen Rendle, Christoph Freudenthaler, Zeno Gantner, and Lars Schmidt-Thieme. 2009.  BPR: Bayesian personalized ranking from implicit feedback. In _UAI_. 452–461.
  * Rosset et al. (2003) Saharon Rosset, Einat Neumann, Uri Eick, and Nurit Vatnik. 2003\.  Customer lifetime value models for decision support.  _Data mining and knowledge discovery_ 7 (2003), 321–339.
  * Sabbeh (2018) Sahar F Sabbeh. 2018\.  Machine-learning techniques for customer retention: A comparative study.  _IJACSA_ 9, 2 (2018).
  * Schmittlein et al. (1987) David C Schmittlein, Donald G Morrison, and Richard Colombo. 1987.  Counting your customers: Who-are they and what will they do next?  _Management science_ 33, 1 (1987), 1–24.
  * Singh et al. (2009) Siddharth S Singh, Sharad Borle, and Dipak C Jain. 2009\.  A generalized framework for estimating customer lifetime value when customer lifetimes are not observed.  _Qme_ 7 (2009), 181–205.
  * Tekin et al. (2022) Ahmet Tezcan Tekin, Tolga Kaya, and Ferhan Cebi. 2022.  Customer lifetime value prediction for gaming industry: fuzzy clustering based approach.  _Journal of Intelligent & Fuzzy Systems_ 42, 1 (2022), 87–96.
  * Tkachenko (2015) Yegor Tkachenko. 2015\.  Autonomous CRM control via CLV approximation with deep reinforcement learning in discrete and continuous action space.  _arXiv preprint arXiv:1504.01840_ (2015).
  * Vanderveld et al. (2016) Ali Vanderveld, Addhyan Pandey, Angela Han, and Rajesh Parekh. 2016.  An Engagement-Based Customer Lifetime Value System for E-commerce. In _KDD_. ACM, 293–302.
  * Venkatesan and Kumar (2004) Rajkumar Venkatesan and Vita Kumar. 2004.  A customer lifetime value framework for customer selection and resource allocation strategy.  _Journal of marketing_ 68, 4 (2004), 106–125.
  * Wang et al. (2022) Fangye Wang, Yingxu Wang, Dongsheng Li, Hansu Gu, Tun Lu, Peng Zhang, and Ning Gu. 2022.  CL4CTR: A Contrastive Learning Framework for CTR Prediction.  _arXiv preprint arXiv:2212.00522_ (2022).
  * Wang et al. (2017) Ruoxi Wang, Bin Fu, Gang Fu, and Mingliang Wang. 2017\.  Deep & cross network for ad click predictions.  In _ADKDD_. 1–7.
  * Wang et al. (2021) Ruoxi Wang, Rakesh Shivanna, Derek Cheng, Sagar Jain, Dong Lin, Lichan Hong, and Ed Chi. 2021.  Dcn v2: Improved deep & cross network and practical lessons for web-scale learning to rank systems. In _WWW_. 1785–1797.
  * Wang et al. (2019) Xiaojing Wang, Tianqi Liu, and Jingang Miao. 2019.  A deep probabilistic model for customer lifetime value prediction.  _arXiv preprint arXiv:1912.07753_ (2019).
  * Wei et al. (2021) Yinwei Wei, Xiang Wang, Qi Li, Liqiang Nie, Yan Li, Xuanping Li, and Tat-Seng Chua. 2021\.  Contrastive learning for cold-start recommendation. In _MM_. 5382–5390.
  * Win and Bo (2020) Than Than Win and Khin Sundee Bo. 2020.  Predicting customer class using customer lifetime value with random forest algorithm. In _ICAIT_. IEEE, 236–241.
  * Wu et al. (2019) Chuhan Wu, Fangzhao Wu, Mingxiao An, Jianqiang Huang, Yongfeng Huang, and Xing Xie. 2019\.  NPA: neural news recommendation with personalized attention. In _KDD_. 2576–2584.
  * Xie et al. (2022) Xu Xie, Fei Sun, Zhaoyang Liu, Shiwen Wu, Jinyang Gao, Jiandong Zhang, Bolin Ding, and Bin Cui. 2022\.  Contrastive learning for sequential recommendation. In _ICDE_. IEEE, 1259–1273.
  * Xing et al. (2021) Mingzhe Xing, Shuqing Bian, Wayne Xin Zhao, Zhen Xiao, Xinji Luo, Cunxiang Yin, Jing Cai, and Yancheng He. 2021\.  Learning Reliable User Representations from Volatile and Sparse Data to Accurately Predict Customer Lifetime Value. In _KDD_. 3806–3816.
  * Yang et al. (2022) Jie Yang, Yilin Li, and Deddy Jobson. 2022.  Personalized Promotion Decision Making Based on Direct and Enduring Effect Predictions.  _arXiv preprint arXiv:2207.14798_ (2022).
  * Yi et al. (2022) Zixuan Yi, Xi Wang, Iadh Ounis, and Craig Macdonald. 2022\.  Multi-modal graph contrastive learning for micro-video recommendation. In _SIGIR_. 1807–1811.
  * Yu et al. (2022) Junliang Yu, Hongzhi Yin, Xin Xia, Tong Chen, Lizhen Cui, and Quoc Viet Hung Nguyen. 2022\.  Are graph augmentations necessary? simple graph contrastive learning for recommendation. In _SIGIR_. 1294–1303.
  * Zhang et al. (2020) Wei Zhang, Xuemei Leng, and Siyu Liu. 2020.  Research on mobile impulse purchase intention in the perspective of system users during COVID-19.  _Personal and Ubiquitous Computing_ (2020), 1–9.
  * Zhao et al. (2023) Shiwei Zhao, Runze Wu, Jianrong Tao, Manhu Qu, Minghao Zhao, Changjie Fan, and Hongke Zhao. 2023.  perCLTV: A general system for personalized customer lifetime value prediction in online games.  _TOIS_ 41, 1 (2023), 1–29.


# Far infrared to terahertz widely tunable narrow linewidth light source via
surface-emitting periodically poled thin film lithium niobate waveguides

Valerie Yoshioka  1 Jicheng Jin  1 Qiang Lin  2 and Bo Zhen 1  1 Department of
Physics and Astronomy, University of Pennsylvania, 209 S 33rd Street,
Philadelphia, PA 19104, USA  2 Department of Electrical and Computer
Engineering, University of Rochester, 500 Wilson Boulevard, Rochester, NY
14627, USA  *<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Generating widely tunable, continuous wave light at long wavelengths via
difference frequency generation (DFG) remains challenging due to high
absorption and dispersion. The relatively new platform of thin film lithium
niobate enables high-confinement nonlinear waveguides, which could improve
efficiency. We simulated DFG in thin film lithium niobate waveguides that are
periodically poled for surface emission at 30 THz. Maximum efficiency for a 1
cm device is $9.16\times 10^{-6}$ W-1 assuming $d_{33}=30$ pm/V. The tuning
range within 50$\%$ of efficiency at 30 THz is as wide as 25 THz, from 25 THz
to 50 THz.

##  1 Introduction

Long wavelength generation is necessary in many different applications,
ranging from bio-sensing [1] to spectroscopy [2] to material detection [3].
While many applications can be performed with long wavelength pulses,
narrowband continuous wave (CW) light tends to require simpler equipment [4]
and can be more easily used for real-time imaging at room temperature [5].
However, efficiently producing widely tunable, narrowband CW light in the far-
infrared (FIR) or terahertz (THz) regime remains difficult due to high
material absorption and rapidly fluctuating dispersion.

Difference frequency generation (DFG) uses the second-order optical
nonlinearity to mix two optical modes that are close in frequency and produce
light at the difference frequency, which can be much longer in wavelength.
While DFG can produce a wide range of wavelengths, phononic resonances in
nonlinear materials introduce high absorption at long wavelengths, often
limiting efficient DFG to short propagation lengths. Various methods to
mitigate this effect have been proposed and can be separated into two
categories: collinear and noncollinear.

In the first category, collinear DFG allows both input optical signals and the
DFG signal to travel together along the nonlinear waveguide. To reduce the
loss felt by the DFG signal, an adjacent waveguide with lower absorption can
be used to guide most of the DFG light [6]. While the long wavelength DFG
signal experiences less loss, it still interacts with the nonlinear material,
allowing the process to continue. However, some portion of the DFG light must
continue to interact with the optical pump signals in the lossy, nonlinear
material for collinear DFG to occur. As such, it experiences exponential loss
as it propagates. While this scheme is advantageous for wavelengths with
reasonable loss in the nonlinear material, it becomes prohibitively
inefficient with higher losses.

In comparison, noncollinear DFG schemes encourage surface emission of long
wavelength light to reduce interaction with the lossy material and thus
prevent absorption. One of the most common noncollinear techniques is
Cherenkov phase matching, which depends on material refractive index. The
Cherenkov condition can be summarized by $\phi_{C}=\arccos\left(\frac{n_{\rm
opt}}{n_{\rm FIR}}\right)$, where $\phi_{C}$ is the emission angle of the DFG
light from the sample, $n_{\rm opt}$ is the refractive index at optical
wavelengths, and $n_{\rm FIR}$ is the refractive index at the longer DFG
wavelength [7]. While this method works effectively in the low THz regime
($<8$ THz) in lithium niobate (LN) due to minimal dispersion, it becomes more
difficult in regions where the refractive index is highly dispersive.

Another method of surface emission interfaces metal antennas with integrated
LN waveguides in order to couple out THz light generated from optical
rectification of pulsed light [8]. While antenna design parameters provide
significant control over the emitted THz waveform, efficiency is limited since
DFG occurs over a short distance, around 200 $\mu$m at most per antenna. As
such, implementing a similar scheme with CW DFG would be inefficient and
limited by these short scales.

An alternate approach uses periodic poling to phase match DFG light to
radiation modes. Instead of being guided in a waveguide, the DFG light
immediately radiates out normal to the waveguide. Analytical derivations of
this behavior [9], along with demonstrations of its performance using PPLN
[10, 11], have been implemented. However, these approaches were limited by
low-confinement, bulk PPLN waveguides [11, 12]. In recent years, thin film
lithium niobate (TFLN) platforms have become commercially available, resulting
in more efficient frequency conversion due to higher confinement in integrated
waveguides [13, 14].

Here, we focus on the potential for using periodically poled TFLN waveguides
as opposed to bulk or diffused waveguides to enhance efficiency and compare
performance with more recent developments in long wavelength generation via
integrated photonic devices. To our knowledge, no integrated TFLN waveguide
has been explored for continuous wave surface emission via periodic poling. In
this work, we demonstrate a design for a surface-emitting, periodically poled
TFLN device that emits widely tunable, narrowband FIR CW light at room
temperature.

##  2 Analytical Model

In a waveguided DFG process, two overlapping optical modes interact in a
nonlinear material, generating a nonlinear polarization that acts as a source
for the desired DFG mode. The frequencies of the optical modes, $f_{1}$ and
$f_{2}$, are chosen such that their difference is the desired output
frequency, $f_{3}=f_{1}-f_{2}$. For x-cut LN waveguides with TE-like optical
modes, the electric field can be written as
$E_{i}=A_{i}e^{i(k_{i}y-\omega_{i}t)}\hat{z}$, where $A_{i}$ is the field
amplitude, $k_{i}=2\pi n_{i}/\lambda_{i}$ is momentum dependent on effective
mode index $n_{i}$ at wavelength $\lambda_{i}$, and $\omega_{i}=2\pi f_{i}$ is
angular frequency. The fields travel along the crystal’s $\hat{y}$ axis and
are directed along its $\hat{z}$ axis. As a result, only the
$\chi^{(2)}_{zzz}=2d_{33}$ term contributes to the nonlinear polarization.
Nonlinear polarization can thus be written as:

$\displaystyle P_{\rm NL}$
$\displaystyle=2\epsilon_{0}\chi^{(2)}_{zzz}E_{1}E_{2}^{*}$ (1) $\displaystyle
P_{\rm NL}$
$\displaystyle=4\epsilon_{0}d_{33}A_{1}A_{2}^{*}e^{i(k_{1}-k_{2})y-i(\omega_{1}-\omega_{2})t}\hat{z}$
(2)

A factor of 2 accounts for both permutations of $E_{1}$ and $E_{2}$. $P_{\rm
NL}$ acts as a source for $E_{3}$ via the wave equation, resulting in the
vector differential equation:

$\displaystyle(\Delta+k_{3}^{2}n_{3}^{2})E_{3}$
$\displaystyle=-k_{3}^{2}\frac{1}{\epsilon_{0}}P_{\rm NL}$ (3)

where $\Delta$ is the vector Laplacian and $n_{3}$ is the refractive index of
LN at $\lambda_{3}$.

Figure 1: DFG with periodic poling for collinear (a-c) and noncollinear (d-f)
approaches. a) Phase matching in the collinear scheme; $k_{1}$ and $k_{2}$ are
the momentum for input optical modes. Poling period $a$ is chosen such that
$k_{1}-k_{2}-k_{3}-\frac{2\pi}{a}=0$. As a result, the $k_{3,x}$ component is
set to 0, and $k_{3}$ is fully directed along the waveguide. b) Nonlinear
polarization $E_{1}\times E_{2}^{*}$ without poling (unshaded) and nonlinear
polarization $P_{\rm NL}$ after poling with collinear period $a$ (shaded). The
generated output field $E_{3,z}$ is a guided mode. c) Power efficiency
$\Gamma$ scales quadratically with length $L$ but decays exponentially with
loss $\alpha$. d) Phase matching in the noncollinear scheme. Poling period $a$
is chosen such that $k_{1}-k_{2}-\frac{2\pi}{a}=0$, so the $k_{3,y}$ component
along the waveguide is set to 0, and $k_{3}$ is fully directed along the
$\hat{x}$ direction. e) Nonlinear polarization $E_{1}\times E_{2}^{*}$ without
poling (unshaded) is the same as the collinear case. Below, nonlinear
polarization $P_{\rm NL}$ after poling (shaded) is shown with the much shorter
noncollinear period $a$ to match the noncollinear condition. $E_{3,z}$ is a
radiating mode that emits normally towards both top and bottom of the
waveguide. c) Efficiency $\Gamma$ scales linearly with length $L$, but there
is no exponential decay with increased length due to loss.

For long enough propagation lengths, output power begins to oscillate rather
than growing due to phase mismatch, $k_{1}-k_{2}-k_{3}\neq 0$. Quasi-phase
matching through ferroelectric periodic poling is often used to counteract
this effect. By flipping the material domain and thus the sign of the
nonlinear term periodically, phase mismatch can be mitigated and the output
power will continue to grow with longer propagation length. This periodic flip
can be described by a Fourier series [15], replacing the constant $d_{33}$ by:

$\displaystyle d(y)$ $\displaystyle=d_{33}\sum_{m}G_{m}e^{ik_{m}y}$ (4)
$\displaystyle G_{m}$
$\displaystyle=\frac{2}{m\pi}\sin\left(\frac{m\pi}{2}\right)$ (5)
$\displaystyle k_{m}$ $\displaystyle=\frac{2\pi m}{a}$ (6)

where $m$ is the index in the Fourier series and $a$ is the poling period. The
largest contributing component occurs for $m=1$ and is most frequently used
for phase matching. With the $m=1$ poling term, the phase gains an additional
term, $\frac{2\pi}{a}$. Using the first and largest poling term, the nonlinear
coefficient can be written as $d(y)=\frac{2}{\pi}d_{33}e^{i\frac{2\pi}{a}y}$.
Substituting in this nonlinear term adjusts the polarization to become:

$\displaystyle P_{\rm NL}$
$\displaystyle=\frac{8}{\pi}\epsilon_{0}d_{33}A_{1}A_{2}^{*}e^{i(k_{1}-k_{2}-\frac{2\pi}{a})y-i(\omega_{1}-\omega_{2})t}\hat{z}$
(7)

Depending on the choice of poling period $a$, the DFG process can be phase
matched such that the majority of output DFG light remains in the waveguide as
a collinear guided mode or emits via a noncollinear radiated mode. In the
collinear case, all $k_{n}$ are directed along $\hat{y}$, and $a$ is chosen
such that $k_{1}-k_{2}-k_{3}-\frac{2\pi}{a}=0$ (Fig. 1a). The resulting
nonlinear polarization gives rise to a guided wave at $\omega_{3}$ (Fig. 1b).
The efficiency scales quadratically with length, but is also exponentially
affected by loss (Fig. 1c).

In comparison, our surface-emitting poling period is chosen to ensure DFG
light immediately leaves the waveguide. In this case, $k_{3,y}$ must vanish
such that $k_{3}$ is directed perpendicular to the waveguide. Thus, we choose
$a$ such that $k_{1}-k_{2}-\frac{2\pi}{a}=0$ exactly (Fig. 1d). In doing so,
the phase-matched process only occurs for modes in which $k_{3,y}=0$. Since
$|k_{3}|$ is defined by the DFG wavelength and material index to be nonzero,
and the electric field $E_{3}$ is polarized in-plane ($\hat{z}$), the DFG
light must radiate normal to the surface of the waveguide (Fig. 1e).
Efficiency scales linearly with length but does not experience exponential
loss with longer propagation lengths since the field quickly radiates away
from the lossy material (Fig. 1f).

Different assumptions are made in solving the wave equation (Eqn. 3) depending
on whether the output is a guided mode or radiated light. In the collinear
case, guided modes can be treated like plane waves with an additional term
that accounts for effective spot area size where all three modes overlap in
the nonlinear material [6]. In the noncollinear case, the output waves radiate
in a cylindrical pattern out from the waveguide, so plane wave assumptions
break down. The analytical model for noncollinear phase matching in a bulk
waveguide has been derived using vector Green functions [9], but our
integrated waveguide model includes a few distinctions.

The first is the inclusion of mode overlap. In a bulk waveguide, the
difference between optical mode sizes is negligible. However, modes in
integrated thin film LN waveguides are more sensitive to changes in
wavelength. Furthermore, while most optical energy is concentrated in the
nonlinear waveguide, some portion extends to the cladding and doesn’t
contribute to DFG; this is especially true for thin, narrow waveguides. When
defining overlap for noncollinear DFG, we only account for overlap between the
two optical modes in the nonlinear waveguide as the DFG field immediately
radiates out of the waveguide. Overlap $\eta$ can thus be written as [12]:

$\displaystyle\eta=\frac{\left(\int_{-w/2}^{w/2}dz\int_{-t/2}^{t/2}E_{1}E_{2}dx\right)^{2}}{\left(\int_{-\infty}^{\infty}dz\int_{-\infty}^{\infty}E_{1}^{2}dx\right)\left(\int_{-\infty}^{\infty}dz\int_{-\infty}^{\infty}E_{2}^{2}dx\right)}$
(8)

where $w$ is the waveguide width and $t$ is the waveguide thickness.

The second difference in our model is simulated collection efficiency. In
comparing our simulation values with the expected efficiency from analytical
models, simulated efficiency is affected by the definition of a measurement
flux plane (Fig. 2b), through which we measure power outflow. In this way, we
include the realistic effects of limited numerical aperture in collection
optics, Fresnel coefficients that depend on dispersion, and losses towards the
substrate. Note that this factor also changes with waveguide parameters since
the energy distribution in the solid angle of the radiated light depends on
waveguide parameters [9].

With these changes, we compare our simulations to the analytical model:

$\displaystyle\Gamma$ $\displaystyle=\frac{P_{3}}{P_{1}P_{2}}=ACL\eta$ (9)
$\displaystyle A$ $\displaystyle=\frac{64\pi
d_{33}^{2}}{\epsilon_{0}cn_{1}n_{2}\lambda_{3}^{3}}$ (10)

where $\Gamma$ is power efficiency, $P_{i}$ is the power at frequency $f_{i}$,
$C$ is the collection efficiency factor obtained by comparing simulated output
power to our analytical model, $L$ is propagation length, $n_{i}$ are
effective mode indices for respective optical modes, $d_{33}$ is the nonlinear
coefficient, and $\eta$ is the mode overlap between the two input optical
modes.

Besides length scaling, noncollinear efficiency differs from collinear
efficiency in a few other critical ways. Noncollinear efficiency scales
cubically with frequency, as opposed to quadratic scaling for collinear
efficiency. As a result, noncollinear efficiency tends to be higher at higher
frequencies, and operation bandwidth is asymmetric. Furthermore, in the
collinear method, effective spot area is important since higher intensity
results in higher conversion efficiency [6]. In comparison, noncollinear
efficiency only depends on normalized optical mode overlap; decreasing mode
area does not necessarily improve efficiency.

Using integrated TFLN waveguides rather than bulk waveguides for noncollinear
DFG still has advantages. Tightly confined optical modes in TFLN waveguides
ensure most electric field energy is present in the nonlinear material for
DFG. Thinner devices also minimize the distance DFG light travels in the
nonlinear material before radiating, reducing loss. Aside from potential
improvements in efficiency, device size is significantly reduced since TFLN
waveguides exhibit better mode confinement.

##  3 Simulation Parameters

Figure 2: a) Diagram of thin film PPLN device. Waveguide parameters include LN
thickness $t$, etch depth $h$, and width $w$. Periodic poling with period $a$
is shown with shading; the direction of the $d_{33}$ component is shown with
black arrows. Optical input light signals $f_{1}$ and $f_{2}$ are guided into
the waveguide. Output light at frequency $f_{3}$ radiates as a cylindrical
wave from the waveguide. b) $E_{3,z}$ radiates from the waveguide as a
cylindrical wave; the labeled flux plane is used to define collected output
power. c) $|E_{1}|$ of TE00 mode for $\lambda_{1}$ = 1.4 $\mu$m with effective
index $n_{1}=2.0229$. Waveguide parameters are $t$ = 800 nm, $h$ = 320 nm, and
$w$ = 2 $\mu$m; sidewall angle is 60°. d) $|E_{2}|$ of TE00 mode for
$\lambda_{2}$ = 1.63 $\mu$m with effective index $n_{2}=1.9843$ for the same
waveguide. e) Waveguide dispersion for guided fundamental TE00 mode below the
light line; $\lambda_{1}$ = 1.4 $\mu$m ($f_{1}$ = 214.3 THz) and $\lambda_{2}$
= 1.63 $\mu$m ($f_{2}$ = 184 THz) are labeled.

Our structure (Fig. 2a) is comprised of a shallow etched TFLN waveguide on top
of SiO2 substrate; the top cladding is air. To demonstrate the capability of
this technique, we simulated emission at $\lambda_{3}$ = 10 $\mu$m ($f_{3}$ =
30 THz), allowing us to fully capture collection efficiency through 3D
simulations with reasonable simulation size and mesh. However, 2D simulations
demonstrate that this technique works with even longer wavelengths into the
THz regime. For waveguide parameters, we used total LN thickness $t$ = 800 nm,
etch depth $h$ = 320 nm, and top width $w$ = 2 $\mu$m with a sidewall angle of
60°. LN waveguides with these dimensions have already been fabricated and
poled, proving their experimental feasibility [16]. Poling period $a$ is set
such that the DFG light emits perpendicular to the top surface unless
otherwise stated.

In this structure, input wavelengths $\lambda_{1}$ = 1.4 $\mu$m and
$\lambda_{2}$ = 1.63 $\mu$m have fundamental TE00 modes with effective
refractive indices of $n_{1}=2.0229$ and $n_{2}=1.9843$ respectively (Fig. 2c
and 2d). In crystal coordinates, TE modes exhibit an electric field pointing
in the $\hat{z}$ direction, while the light propagates along the waveguide in
the $\hat{y}$ direction. Waveguide dispersion shows that both of these TE00
modes are guided as they lie below the light line (Fig. 2e).

We used 30 pm/V as the $d_{33}$ value for our DFG simulations. While there is
evidence that $d_{33}$ can be much larger around 181-195 pm/V in the THz
regime [17, 6], absolute $\chi^{(2)}$ measurements involving frequency mixing
between vastly different wavelengths are limited. Furthermore, the
$\chi^{(2)}$ values are expected to dramatically fluctuate in the the
Reststrahlen regime due to strong dispersion. For LN in the region between
1-100 THz, $\chi^{(2)}_{zzz}$ is expected to vary between values as high as
2000 pm/V to values as low as 2 pm/V [18]. However, this dramatic fluctuation
has not been rigorously verified experimentally. As such, we set $d_{33}$ to
be 30 pm/V as a characteristic value. Note that any changes in this value due
to dispersion will cause efficiency to scale proportionally to $d_{33}^{2}$.
If we used $d_{33}$ = 195 pm/V for example, our efficiency would scale up by a
factor of 42. In general, additional experimental data on $\chi^{(2)}$ values
in the THz regime would be helpful to confirm nonlinear strength at particular
wavelengths.

We used COMSOL to simulate the surface emitting DFG process in periodically
poled TFLN waveguides. Input modes were simulated and injected into the
waveguide. These optical fields were used to define the nonlinear
polarization, which acted as the source for the DFG light at frequency
$f_{3}$. We then simulated the resulting DFG radiation as a cylindrical wave
originating from the waveguide (Fig. 2b). To determine expected efficiency, we
evaluated the surface integral for power emitted through the flux plane at the
top surface. The collection efficiency factor therefore includes losses due to
emission towards the substrate and sides, which better reflects the expected
experimental efficiency.

##  4 Results

In this scheme, efficiency scales linearly with power (Fig. 3). If the
structure is too short, around $L<10a$, the efficiency is lower than expected
due to less uniform power distribution from edge effects. However, as length
is increased and edge effects diminish, efficiency trends towards linear
increase with propagation length as expected by the analytical model.
Reasonable waveguide lengths are much greater than this transition point, so
experimental devices will exhibit a linear increase in efficiency with length.
For a 1 cm device with 30 THz output, we expect efficiency around $9.16\times
10^{-6}$ W-1.

Figure 3: Efficiency v. length. In general, efficiency scales linearly with
length. The deviation at much smaller lengths is due to edge effects
dominating the signal. Figure 4: DFG emission with varying poling period. a)
Normalized efficiency v. poling period; efficiency is normalized by maximum
efficiency, which occurs at a poling period of 4.4 $\mu$m (red circle). b)
Emission angle v. poling period; the poling period of 4.4 $\mu$m (red circle)
has an emission angle around 0°. c) Emission directions for different poling
periods $a$ = 4.2 $\mu$m, 4.4 $\mu$m, and 4.6 $\mu$m. Emission angle changes
as poling period shifts. Note that emission towards the substrate dies out
quickly due to higher loss in oxide at 30 THz.

When poling period is adjusted, radiated output light emits at an angle that
deviates from fully perpendicular (Fig. 4). This emission angle $\theta$ can
be controlled by choosing $a$ such that
$k_{1}-k_{2}-\frac{2\pi}{a}=k_{3,y}=k_{3}\sin\theta$. $\theta$ can be positive
or negative depending on emission direction but cannot exceed 90° in either
direction, limiting the range for poling periods. For our device, poling
period must be within the range 2.8 $\mu$m < $a$ < 10.6 $\mu$m to allow
emission at 30 THz.

Efficiency is highest for normal emission ($\theta$ = 0°, $a$ = 4.4 $\mu$m)
since our efficiency is defined by evaluating output power through the top
surface. Normal emission also reduces loss due to Fresnel reflection
coefficients. A steeper angle can result in internal reflection, reducing
outgoing light. Fresnel coefficients change according to material dispersion,
so critical angles for total internal reflection will depend on the desired
central frequency and tuning range.

Figure 5: Normalized efficiency v. output frequency for a fixed waveguide
structure with $a$ = 4.4 $\mu$m. Efficiency is normalized relative to the
efficiency at $f_{3}$ = 30 THz (red circle). Higher frequencies are more
efficient due to cubic dependence of efficiency on frequency. The large peak
around 35 THz is due to a phonon resonance in the SiO2 substrate, which
effectively eliminates substrate loss.

To determine bandwidth for a single device, we fixed the poling period to be
4.4 $\mu$m and swept output frequency (Fig. 5). Material dispersion of both LN
and the SiO2 substrate was included, resulting in a large peak around 35 THz.
At this frequency, phonon resonances in SiO2 affect its refractive index such
that imaginary index is larger than real index. The substrate therefore acts
like a metal near this frequency, eliminating substrate loss and greatly
increasing efficiency. In general, higher frequencies with similar loss tend
to be more efficient due to the cubic dependence of efficiency on output
frequency $f_{3}=c/\lambda_{3}$. The bandwidth around 30 THz, here defined as
the full width at half of the efficiency at 30 THz, ranges from about 25 THz
to 50 THz. While asymmetrical, the overall bandwidth spans 25 THz, allowing
wide tuning of output wavelengths for different applications.

##  5 Comparison to Other Sources

In collinear DFG, output power scales quadratically with propagation length.
However, collinear DFG is exponentially affected by absorption. Even though
there are clever approaches to reduce effective loss that the DFG mode
experiences, such as guiding it in an alternate material with lower loss [6],
some portion of the DFG light still needs to travel in the lossy, nonlinear
material to benefit from quadratic scaling of collinear scheme, experiencing
exponential loss. After a certain length, the exponential absorption is
stronger than the quadratic growth with length, and efficiency decreases with
longer devices (Fig. 6); this is especially problematic at frequencies above 5
THz where imaginary refractive index $\kappa$ spikes due to phonon resonances
in LN. When $\kappa$ is large, the length at which maximum efficiency occurs
is too short to efficiently generate long wavelength light.

Figure 6: Efficiency v. propagation length in both collinear (solid line) and
noncollinear (dotted line) schemes for different losses set by imaginary
refractive index $\kappa$ (color-coded). Output frequency $f_{3}$ is 30 THz,
and $d_{33}$ is 30 pm/V for both schemes. Differences in noncollinear
efficiency at $\kappa=0,0.01,0.1$ are negligible on this scale. For collinear
efficiency, phase mismatch affects efficiency through a sinc squared term,
resulting in oscillations starting around 5 mm for $\kappa=0$. For losses of
$\kappa=0.01$ and above, absorption begins to dominate, washing out effects of
phase mismatch.

In comparison, our noncollinear scheme exhibits efficiency that continues to
scale linearly with length regardless of loss since DFG light is immediately
emitted and experiences a constant level of absorption. The efficiency per
length slightly decreases with higher loss, but propagation length is no
longer limited. For long enough propagation lengths, noncollinear devices
become more efficient than collinear devices; this crossover length depends on
the absorption of the material at the desired wavelength. At 30 THz with
$\kappa=0.0431$ for LN [18], the noncollinear case is more efficient after
propagation lengths above 1 mm (Fig. 6), which is reasonable for integrated
optics. For $\kappa=1$, a noncollinear device is more efficient than a
collinear one regardless of length.

Furthermore, collinear schemes are subject to strict phase matching demands,
and after a certain length, efficiency is saturated and begins to oscillate.
In comparison, noncollinear phase matching requirements are less stringent
since a slight deviation in phase matching simply matches the output light to
an angled radiation mode, allowing DFG to continue. Efficiency will continue
to increase with propagation length, until absorption or depletion of input
wavelengths start to become significant.

In comparison to Cherenkov surface emission schemes, our design is less
affected by strong dispersion. In bulk devices, Cherenkov emission is not
possible between 10-15 THz or >20 THz since LN refractive index for THz light
dips below that for optical light. In waveguide devices, engineering waveguide
parameters and using other adjacent materials can tune the effective index for
THz light, allowing for Cherenkov emission at otherwise impossible frequencies
[18]. However, Cherenkov surface emission for continuous wave DFG has only
been demonstrated in LN ion-planted slab and channel waveguides with low index
contrast [7], resulting in relatively weak confinement. Furthermore, many
devices that rely on Cherenkov emission require silicon prisms on top of the
nonlinear material to improve emission towards the top surface since Fresnel
coefficients cause internal reflection at the Cherenkov angle [18, 7]. Since
emission angle in our approach relies solely on the choice of poling period,
we have a greater degree of freedom in our design to control emission angle in
order to maximize efficiency for a central wavelength.

A fair comparison between our simulated efficiency and previous work on
surface emission based on periodic poling is difficult since these works used
different output wavelengths, nonlinear strength, and underlying assumptions
[9, 12]. Since many of these parameters are also wavelength dependent, it
becomes even more difficult to separate differences in models from how well
efficiency is optimized in the waveguide structure. That said, our work uses
simulations rather than relying solely on analytical calculations; it thus
includes many contributing effects, ranging from material dispersion, mode
overlap, limited collection angle, Fresnel coefficients, and any other
physical effects that are difficult to express analytically.

##  6 Conclusion

We have simulated surface emission of 10 $\mu$m CW light using periodically
poled TFLN waveguides and demonstrated their potential for efficient
generation of narrowband, widely tunable CW light at long wavelengths. While
material dispersion still affects device performance at different frequencies,
efficiency scales linearly with length and is not exponentially affected by
loss. Output power can continually be improved by increasing device length
until other effects like pump depletion arise. Furthermore, our device remains
reasonably efficient in generating CW light rather than requiring strong
pulses, enabling production of narrowband, long wavelength light for
applications that require a more precise wavelength or CW light. Frequency
tuning also continuously adjusts emission angle, potentially enabling beam
steering. In the future, we will consider phase modulating an array of
adjacent waveguides to beam steer long wavelength light, allowing even better
control over emission angle and tunability.

Funding This work was partly supported by the U.S. Office of Naval Research
(ONR) through grant N00014-21-1-2703, the Army Research Office under award
contract W911NF-19-1-0087, and DARPA under agreement HR00112220013.

## References

  * [1] P. Lucas, Z. Yang, M. K. Fah, T. Luo, S. Jiang, C. Boussard-Pledel, M.-L. Anne, and B. Bureau, “Telluride glasses for far infrared photonic applications,” Optical Materials Express 3, 1049 (2013).
  * [2] A. Wulf, K. Fumino, R. Ludwig, and P. F. Taday, “Combined THz, FIR and Raman Spectroscopy Studies of Imidazolium-Based Ionic Liquids Covering the Frequency Range 2-300 cm-1,” ChemPhysChem 11, 349–353 (2010).
  * [3] W. H. Fan, A. Burnett, P. C. Upadhya, J. Cunningham, E. H. Linfield, and A. G. Davies, “Far-Infrared Spectroscopic Characterization of Explosives for Security Applications Using Broadband Terahertz Time-Domain Spectroscopy,” Appl. Spectrosc. 61, 638–643 (2007).
  * [4] N. Karpowicz, H. Zhong, J. Xu, K. I. Lin, J. S. Hwang, and X. C. Zhang, “Comparison between pulsed terahertz time-domain imaging and continuous wave terahertz imaging,” Semiconductor Science and Technology 20 (2005).
  * [5] Y. Zhang, C. Wang, B. Huai, S. Wang, Y. Zhang, D. Wang, L. Rong, and Y. Zheng, “Continuous-Wave THz Imaging for Biomedical Samples,” Applied Sciences (Switzerland) 11, 1–26 (2021).
  * [6] J. Yang and C. Wang, “Efficient terahertz generation scheme in a thin-film lithium niobate-silicon hybrid platform,” Optics Express 29, 16477 (2021).
  * [7] M. De Regis, L. Consolino, S. Bartalini, and P. De Natale, “Waveguided Approach for Difference Frequency Generation of Broadly-Tunable Continuous-Wave Terahertz Radiation,” Applied Sciences 8, 2374 (2018).
  * [8] A. Herter, A. Shams-Ansari, F. F. Settembrini, H. K. Warner, J. Faist, M. Lončar, and I. C. Benea-Chelmus, “Terahertz waveform synthesis in integrated thin-film lithium niobate platform,” Nature Communications 14 (2023).
  * [9] T. Suhara, Y. Avetisyan, and H. Ito, “Theoretical analysis of laterally emitting terahertz-wave generation by difference-frequency generation in channel waveguides,” IEEE Journal of Quantum Electronics 39, 166–171 (2003).
  * [10] K. Suizu, Y. Suzuki, Y. Sasaki, H. Ito, and Y. Avetisyan, “Surface-emitted terahertz-wave generation by ridged periodically poled lithium niobate and enhancement by mixing of two terahertz waves,” Optics Letters 31, 957 (2006).
  * [11] J. Hamazaki, Y. Ogawa, T. Kishimoto, S. Hayashi, N. Sekine, and I. Hosako, “Conversion efficiency improvement of terahertz wave generation laterally emitted by a ridge-type periodically poled lithium niobate,” Optics Express 30, 11472 (2022).
  * [12] Y. J. Ding, “Efficient generation of far-infrared radiation from a periodically poled LiNbO3 waveguide based on surface-emitting geometry,” Journal of the Optical Society of America B 28, 977 (2011).
  * [13] G. Chen, N. Li, J. D. Ng, H. L. Lin, Y. Zhou, Y. H. Fu, L. Y. T. Lee, Y. Yu, A. Q. Liu, and A. J. Danner, “Advances in lithium niobate photonics: development status and perspectives,” Advanced Photonics 4, 1–43 (2022).
  * [14] A. Boes, L. Chang, C. Langrock, M. Yu, M. Zhang, Q. Lin, M. Lončar, M. Fejer, J. Bowers, and A. Mitchell, “Lithium niobate photonics: Unlocking the electromagnetic spectrum,” Science (New York, N.Y.) 379, eabj4396 (2023).
  * [15] R. W. Boyd, _Nonlinear Optics_ (Academic Press, 1992), 3rd ed.
  * [16] X. Wang, X. Jiao, B. Wang, Y. Liu, X.-P. Xie, M.-Y. Zheng, Q. Zhang, and J.-W. Pan, “Quantum frequency conversion and single-photon detection with lithium niobate nanophotonic chips,” npj Quantum Information 9, 38 (2023).
  * [17] G. D. Boyd and M. A. Pollack, “Microwave nonlinearities in anisotropic dielectrics and their relation to optical and electro-optical nonlinearities,” Physical Review B 7, 5345–5359 (1973).
  * [18] B. N. Carnio and A. Y. Elezzabi, “Investigation of ultra-broadband terahertz generation from sub-wavelength lithium niobate waveguides excited by few-cycle femtosecond laser pulses,” Optics Express 25, 20573 (2017).