# Analyzing performance of $f$-deformed displaced Fock state in continuous-
variable quantum teleportation

Deepak<EMAIL_ADDRESS>Arpita Chatterjee<EMAIL_ADDRESS>Department of Mathematics
J. C. Bose University of Science and Technology, YMCA
Faridabad 121006, Haryana, India

###### Abstract

Performing non-Gaussian operations, namely photon addition, photon
subtraction, photon-addition-then-subtraction, photon-subtraction-then-
addition can successfully enhance the fidelity of the continuous-variable
quantum teleportation. However, a shortcoming of these non-Gaussian resources
is that they are probabilistic in nature. In this article, we investigate the
success probability of the non-Gaussian resources for optimal performance of
the ideal teleportation protocol. To this end, we first derive the analytical
expression for the two-mode entangled channel having $f$-deformed displaced
Fock state or photon-added displaced Fock state or photon-subtracted displaced
Fock state at one port and vacuum at another port of a symmetric beam-
splitter. The generalized displaced Fock states are obtained by replacing the
conventional bosonic functions by the nonlinear $f$-deformed operators such as
$A=af(a^{\dagger}a)$ and $B=af(a^{\dagger}a)^{-1}$. The Wigner characteristic
functions describing these three aforementioned non-Gaussian states are
determined and utilized to attain the teleportation fidelity for input
coherent and squeezed vacuum states. It is found that the nonlinear
substitution leads to an enhancement in teleportation fidelity beyond the
threshold limit. Moreover, the entangled photon-subtracted displaced Fock
state demonstrates maximum efficiency as a quantum channel for teleporting
single-mode coherent and squeezed states. We provide the squeezing regime
values corresponding to the optimal performance of the non-Gaussian states
considered, which will be of significant interest to the experimental
fraternity. Further, we show that the entangled photon-added displaced Fock
states have larger amount of entanglement but the entangled photon-subtracted
displaced Fock states have stronger Einstein-Podolsky-Rosen correlation.

$f$-deformed displaced Fock state, non-Gaussian operation, continuous-variable
quantum teleportation

###### pacs:

42.50.-p, 42.50.Dv

##  I Introduction

It is well known that Gaussian type entangled resource such as a two-mode
squeezed vacuum state (TMSV) and Gaussian operations implemented by any
combination of linear, quadratic optics, and homodyne detection are generally
employed in continuous-variable quantum information processing tasks [1, 2].
However, it has been realized that both Gaussian operations and entangled
resources have some limitations. For example, quantum speed up for harmonic
oscillators cannot be attained only by Gaussian operations with Gaussian
inputs [3, 4], entanglement distillation from two Gaussian entangled states is
not possible only by Gaussian local operations and classical communication [5,
6, 7]. Moreover, it has been proved recently that using Gaussian operations in
quantum communication protocols cannot protect Gaussian states from Gaussian
errors [8]. Therefore, it becomes essential to look for non-Gaussian resources
and operations which are more efficient in quantum information processing. One
possible method for producing non-Gaussian entangled resources is to perform
photon addition and/or subtractio on a given Gaussian state [1, 2]. Current
advancements in experimental techniques have enabled the performance of photon
addition and/or subtraction on a given input state.

In recent time, Priya et. al. [9] studied the nonclassical properties of a set
of engineered quantum states (both photon added and photon subtracted) using
the lower and higher-order nonclassicality criteria. In some other work, they
have quantified the nonclassicality and non-Gaussianity of photon-added
displaced Fock states and discussed their significance in quantum technology
and the necessity to quantify these resources [1]. An entangled two-mode state
can be originated by employing a 50:50 symmetric beam-splitter (BS) while a
Fock state and a vacuum are injected at the input modes of the BS.
Consequently, inserting a class of generalized Fock states, e.g. displaced
Fock state (DFS), photon-added displaced Fock state (PADFS), and photon-
subtracted displaced Fock state (PSDFS) etc. from one input port and vacuum
from the other of a 50:50 BS can manufacture two-mode entangled DFS (TMEDFS),
two-mode entangled PADFS (TMEPADFS), and two-mode entangled PSDFS (TMEPSDFS)
respectively. Deepak et al. [2] performed a quantitative investigation of the
nonclassical and quantum non-Gaussian characters of the photon-subtracted
displaced Fock stste by using a collection of measures like Wigner logarithmic
negativity, linear entropy potential, skew information-based measure, and
relative entropy of quantum non-Gaussianity. In present work, the performance
of the TMEPADFS and TMEPSDFS as a quantum channel for continuous-variable
quantum teleportation is investigated. In this approach, the input state is
assumed to be either coherent or squeezed, and the success under the
Braunstein and Kimble protocol is evaluated in terms of the average fidelity
for transporting that single-mode optical state. At the same time, we work to
comprehend how various realistic parameters contribute to successful
teleportation of a single-mode input state. The advantage of $f$-deformed
entangled resources in continuous-variable treleportation protocol under ideal
circumstances is also explored.

The idea of quantum teleportation was first proposed by Bennett et al. in the
context of discrete variables [10], and later illustrated experimentally by
Bouwmeester [11] and Popescu [12]. Vaidman introduced the concept of quantum
teleportation in the continuous-variable regime [13]. The actual quantum-
optical protocol for transmitting quadrature amplitudes of a light field from
a sender to a receiver in a distant place was developed by Braunstein and
Kimble [14] and implemented experimentally by Furusawa et al. [15, 16, 17]
soon afterward. In a typical teleportation protocol, two users, Alice and Bob,
shared an entangled resource of two modes $A$ and $B$. The state to be
teleported is a single-mode input state, called $\ket{\text{in}}$ in the
sender’s place. The input mode $\ket{\text{in}}$ and the mode $A$ of the
entangled channel are mixed at a 50:50 beam-splitter yielding the output modes
$\ket{\text{out}}$ and $A^{\prime}$. Then Alice performed a destructive
homodyne measurement on the output modes and classically communicated the
result to Bob. Bob applied a unitary displacement operation on mode $B$ to
receive the desired teleported state. Different alternative descriptions of
the original Braunstein-Kimble protocol have been used in the literature [18,
19, 20]. The quantum information fraternity is particularly interested in
usage of non-Gaussian resources for improving the theoretical description of
CV quantum teleportation and its experimental implementation due to the
escalating application of CV entangled resources in quantum information
processing and quantum computing [21]. Very recently, Deepak et al. [22]
described the ideal and realistic continuous-variable quantum teleportation
protocol in terms of the characteristic functions of the quantum states
involved (input, resource, and teleported states). Here we have extended that
work by employing the generalized two-mode $f$-deformed entangled resources as
quantum channels.

Entanglement and Einstein-Podolsky-Rosen (EPR) correlation are fundamental
phenomena in the context of quantum mechanics. These concepts probe into the
mysterious and paradoxical nature of the quantum world, challenging our
classical understanding of reality. Entanglement refers to a unique quantum
phenomenon where two or more particles become interconnected in such a way
that the state of one particle instantaneously influences the state of
another, irrespective of the physical distance between them. This correlation
exists even when the particles are in remote location, defying our classical
notions of causality and locality. EPR correlation, proposed by Einstein,
Podolsky and Rosen, is an important aspect of entanglement. In their
pioneering research [23], these physicists proposed a hypothetical experiment
to illustrate what they have noticed as a paradox within quantum mechanics.
They have found that two quantum-entangled particles can have perfectly
correlated positions and momenta that means the measurement of one entangled
particle’s property could immediately determine the state of another, that
violates the principle of local realism where nothing can travel faster than
the speed of light. In this paper, the entanglement amount and the EPR
correlation of the non-Gaussian states of interest are investigated. These
states are obtained by performing the nonlinear $f$-deformation and mixing a
photon-added/subtracted displaced Fock state with a vacuum at a 50:50 beam-
splitter. The performance of these non-Gaussian entangled resources in
continuous-variable quantum teleportation of Braunstein and Kimble protocol
for coherent states, and squeezed states is studied.

The paper is structured as follows: the characteristic functions of all
quantum states under consideration are studied in Sec. II. In Sec. III, we
have calculated the fidelity for continuous-variable quantum teleportation
based on the non-Gaussian entangled resources. In Sec. IV, the entanglement
and EPR correlation are obtained. In the last Sec. V, we have summarized the
main results.

##  II $f$-deformed Displaced Fock States

A displaced Fock state (DFS) described as a generalized coherent state, is
defined by $\ket{\psi}=D(\alpha)\ket{n}$ where $D(\alpha)$ is the displacement
operator and $\ket{n}$ is an arbitrary Fock state. In special case ($n=0$),
this state reduces to a displaced vacuum state, also known as coherent state.
If $\alpha=0$, the displaced Fock state represents a number state.

The wave functions for coherent, squeezed, $k$-photon-added and $k$-photon-
subtracted displaced Fock states may be described as

$\displaystyle\ket{\psi_{\textrm{coh}}}$
$\displaystyle=\ket{\alpha_{0}}=D(\alpha_{0})\ket{0}$
$\displaystyle\ket{\psi_{\textrm{squ}}}$
$\displaystyle=\ket{z}=S(z)\ket{0},\,\,z=re^{i\zeta}$ (1)
$\displaystyle\ket{\psi_{+}}$
$\displaystyle=N_{+}a^{{\dagger}k}D(\alpha)\ket{n}$
$\displaystyle\ket{\psi_{-}}$ $\displaystyle=N_{-}a^{k}D(\alpha)\ket{n}$

Here subscripts $\pm$ correspond to photon addition (photon subtraction),
$N^{-2}_{+}=\sum_{p=0}^{k}{k\choose p}^{2}|\alpha|^{2(k-p)}\frac{(n+p)!}{n!}$
and $N^{-2}_{-}=\sum_{p=0}^{k}{k\choose
p}^{2}|\alpha|^{2(k-p)}\frac{n!}{(n-p)!}$ are the normalization constants for
$k$-photon-added displaced Fock state (PADFS) and $k$-photon-subtracted
displaced Fock state (PSDFS), respectively. These states can be expressed in
terms of the Fock state basis as following:

$\ket{\psi_{\pm}}=\sum_{m=0}^{\infty}C_{m}^{\pm}(n,k,\alpha)\ket{m}$ (2)

where $C_{m}^{+}(n,k,\alpha)=N_{+}\bra{m}a^{{\dagger}k}D(\alpha)\ket{n}$ and
$C_{m}^{-}(n,k,\alpha)=N_{-}\bra{m}a^{k}D(\alpha)\ket{n}$. To investigate the
performance of two classes of non-Gaussian entangled resources for
teleportation of input coherent and squeezed states, we require the two-mode
entangled channel $\rho_{AB}$ that is originated by mixing a PADFS (PSDFS) and
a vacuum at a 50:50 beam-splitter (BS). The symmetric beam-splitter performs a
balanced mixing of a Fock state and a vacuum injected at its input ports
through a $SU(2)$ transformation that yields [2]:

$\ket{n}\otimes\ket{0}=\ket{n,0}\xmapsto{\textbf{BS}}\frac{1}{2^{n/2}}\sum_{j=0}^{n}\sqrt{{n\choose
j}}\ket{j,n-j}$ (3)

Thus inserting PADFS (PSDFS) and vacuum at the input modes of the 50:50 beam-
splitter, the post-BS two-mode entangled photon-added displaced Fock state
(TMEPADFS) (two-mode entangled photon-subtracted displaced Fock state
(TMEPSDFS)) $\rho_{AB}=\ket{\phi_{\pm}}\bra{\phi_{\pm}}$ can be obtained as

$\displaystyle\ket{\phi_{\pm}}$ $\displaystyle=\ket{\psi_{\pm}}\otimes\ket{0}$
$\displaystyle\xmapsto{\textbf{BS}}\sum_{m=0}^{\infty}\frac{C_{m}^{\pm}(n,k,\alpha)}{2^{m/2}}\sum_{j=0}^{m}\sqrt{m\choose
j}\ket{j,m-j}$ (4)

In order to introduce nonlinearity to TMEPADFS (TMEPSDFS), the generalized
displacement operator $D^{\prime}_{f}(\alpha)$ may be of use. This can be
accomplished by recalling the first postulate of the nonlinear coherent states
approach by Man’ko et. al. [24]. The postulate affirms that the standard
bosonic operators undergo $f$-deformation by applying an intensity-dependent
function as follows:

$\displaystyle A$ $\displaystyle=af(N)=f(N+1)a,$ $\displaystyle A^{\dagger}$
$\displaystyle=f^{\dagger}(N)a^{\dagger}=a^{\dagger}f^{\dagger}(N+1)$

Here $f(N)$ is a number operator-valued ($N\equiv a^{\dagger}a$) nonlinear
function of the intensity of light. Considering the $f$-deformed (nonlinear)
bosonic operators, the generalized displacement operator $D_{f}(\alpha)$ can
be written as

$D_{f}(\alpha)=\exp(\alpha A^{\dagger}-\alpha^{*}A)$ (5)

For further expansion of this generalized displacement operator, it is
necessary to find a way to disentangle the $f$-deformed bosonic operators
using Baker-Campbell-Hausdorf formula. The method of auxiliary operators
suggests two auxiliary operators, namely $B=a{f(N)}^{-1}$ and
$B^{\dagger}={f^{{\dagger}}(N)}^{-1}a^{\dagger}$ to overcome this difficulty
[25]. This leads to the generators $\\{A,B^{\dagger},B^{\dagger}A,I\\}$ and
$\\{B,A^{\dagger},A^{\dagger}B,I\\}$ that satisfy the Weyl–Heisenberg Lie
algebra as well as the commutation relations $[A,~{}B^{\dagger}A]=A$,
$[B^{\dagger},~{}B^{\dagger}A]=-B^{\dagger}$, and
$[A,~{}B^{\dagger}]=[B,~{}A^{\dagger}]=I$. As a consequence of the auxiliary
operators method, two generalized displacement operators are defined as

$\displaystyle D^{\prime}_{f}$ $\displaystyle=\exp(\alpha
A^{\dagger}-\alpha^{*}B)$ $\displaystyle=\exp(-|\alpha|^{2}/2)\exp(\alpha
A^{\dagger})\exp(-\alpha^{*}B)$ (6) $\displaystyle D^{\prime\prime}_{f}$
$\displaystyle=\exp(\alpha B^{\dagger}-\alpha^{*}A)$
$\displaystyle=\exp(-|\alpha|^{2}/2)\exp(\alpha
B^{\dagger})\exp(-\alpha^{*}A)$ (7)

Considering the nonlinear displacement operator $D^{\prime}_{f}(\alpha)$, the
generalized $f$-deformed two-mode entangled displaced Fock resource states can
be expressed as

$\ket{\varphi_{\pm}}=\sum_{m=0}^{\infty}\frac{\mathbb{C}_{m}^{\pm}(n,k,\alpha)}{2^{m/2}}\sum_{j=0}^{m}\sqrt{m\choose
j}\ket{j,m-j}$ (8)

where
$\mathbb{C}_{m}^{+}(n,k,\alpha)=\mathbb{N}_{+}\bra{m}A^{{\dagger}k}D^{\prime}_{f}(\alpha)\ket{n}$
and
$\mathbb{C}_{m}^{-}(n,k,\alpha)=\mathbb{N}_{-}\bra{m}A^{k}D^{\prime}_{f}(\alpha)\ket{n}$,
$\mathbb{N}_{\pm}$ are the corresponding normalization constants.

In the continuous-variable teloportation following the Braunstein-Kimble (BK)
protocol, the fidelity is described in terms of the characteristic functions
of the quantum states involved (input, resource, and teleported states) [26].
This formalism is especially beneficial while dealing with non-Gaussian states
and resources because it very much simplifies calculation complexities. The
characteristic functions for a coherent state $\ket{\alpha_{0}}$, a squeezed
state $\ket{z}$, and a TMEPADFS (TMEPSDFS) $\ket{\phi_{\pm}}$ can be obtained
as

$\displaystyle\chi_{\textrm{coh}}(\gamma)$
$\displaystyle=\exp(\frac{-|\gamma|^{2}}{2}+\alpha_{0}^{\star}\gamma-\alpha_{0}\gamma^{\star}),$
(9) $\displaystyle\chi_{\textrm{squ}}(\gamma)$
$\displaystyle=\exp(-\frac{2|\gamma|^{2}\cosh{2r}+(\gamma^{2}e^{-\iota\zeta}+\gamma^{\star
2}e^{\iota\zeta})\sinh{2r}}{4})$ (10)

and

$\displaystyle\chi_{12}^{\pm}(\gamma_{1},\gamma_{2})$
$\displaystyle=\exp(-|\gamma_{1}|^{2}/2-|\gamma_{2}|^{2}/2)$
$\displaystyle\times\sum_{m,l=0}^{\infty}\frac{\left(C_{m}^{\pm}(n,k,\alpha)\right)^{*}C_{l}^{\pm}(n,k,\alpha)}{2^{(m+l)/2}}\sqrt{\frac{l!}{m!}}\gamma_{2}^{m-l}$
$\displaystyle\times\sum_{j,k=0}^{m,l}{m\choose
j}\gamma_{1}^{j-k}\gamma_{2}^{k-j}L_{k}^{j-k}(|\gamma_{1}|^{2})L_{l-k}^{m-l+k-j}(|\gamma_{2}|^{2})$
(11)

where $L_{n}^{k}(x)$ is the associated Laguerre polynomial of order $n$.

##  III Teleportation using $f$-deformed two-mode entangled resources

The idea of quantum teleportation was first proposed by Bennett et al. in the
discrete variable regime [10]. Later on, Braunstein and Kimble (BK) introduced
the actual continuous-variable (CV) quantum-optical protocol for teleporting
the quadrature amplitudes of a light field that employed the Wigner functions
of the associated quantum states [14]. For any given state, there is a one-to-
one correspondence between its Wigner function and symmetrically ordered
characteristic function. The BK protocol in terms of the characteristic
functions is particularly suitable while employing non-Gaussian states as it
greatly simplifies the calculation complexities.

In the standard BK protocol of continuous-variables, two users, Alice and Bob
shared a two-mode entangled state of modes $A$ and $B$. A single-mode input
quantum state to be teleported is provided to Alice. The characteristic
functions of the single-mode input state to be teleported and the two-mode
entangled resource state are denoted by $\chi_{\textrm{in}}(\gamma)$ and
$\chi_{\textrm{res}}(\gamma^{*},\gamma)$, respectively. The input state and
mode $A$ of the entangled resource are mixed at a 50:50 beam-splitter. After
that, the sender performed a destructive homodyne measurement at the output
modes of the BS and the results are classically communicated to the receiver.
Based on Alice’s outcome, Bob applied a specific set of quantum gates on mode
$B$ of the entangled pair and displayed his result leading to the teleported
state. By exploiting the Weyl expansion, the BK protocol allows one to write
the characteristic function of the teleported state in a factorized form as

$\chi_{\textrm{out}}(\gamma)=\chi_{\textrm{in}}(\gamma)\chi_{\textrm{res}}(\gamma^{*},\gamma)$
(12)

In order to measure the probability of success in a teleportation protocol, it
is convenient to use the teleportation fidelity $F$. This is a state-dependent
quantity that measures the overlap between the input state
$\rho_{\textrm{in}}$ and the output (teleported) state $\rho_{\textrm{out}}$,
i.e., $F=\textrm{Tr}[\rho_{\textrm{in}}\rho_{\textrm{out}}]$. If the input and
output states are identical, then $F=1$ corresponding that perfect
teleportation has been reached [27]. The fidelity in terms of the associated
characteristic functions reads as [28, 29]

$F=\frac{1}{\pi}\int
d^{2}\gamma\chi_{\textrm{in}}(\gamma)\chi_{\textrm{out}}(-\gamma),$ (13)

It has been found that a maximum fidelity of 1/2 can be achieved when an input
coherent state is teleported without using a shared entangled state.
Therefore, a successful quantum teleportation is observed when the fidelity
value surpasses the classical limit of 1/2. We require an infinitely entangled
resource state to achieve perfect teleportation with unit fidelity. Further,
the fidelity should be greater than 2/3 for secure teleportation, i.e., to
ensure that any copy of the teleported state is not possible by an
eavesdropper [30].

In the following, we will make use of (13) to derive the fidelity for
teleportation of input coherent and squeezed vacuum states and to quantify the
success of the protocol for different classes of non-Gaussian entangled
resources.

###  III.1 Input Coherent State

The fidelity for teleporting an input coherent state via TMEPADFS (TMEPSDFS)
resource states can be evaluated using (13) as

$\displaystyle F({\phi_{\pm}^{\textrm{coh}}})$
$\displaystyle=\sum_{m,l=0}^{\infty}\frac{(C_{m}^{\pm}(n,k,\alpha))^{*}C_{l}^{\pm}(n,k,\alpha)}{2^{m+1}}\sqrt{\frac{l!}{m!}}(-1)^{m-l}$
$\displaystyle\times\sum_{j=0}^{m}{m\choose
j}\sum_{r,s=0}^{j+\frac{l-m}{2},-j+\frac{l+m}{2}}{j\choose j+\frac{l-m}{2}-r}$
(14)
$\displaystyle\times{m-j\choose-j+\frac{l+m}{2}-s}\frac{(-1)^{r+s}}{2^{r+s}r!s!}\left(\frac{m-l}{2}+r+s\right)!$

Figure 1: (Color online) Variation of fidelity
$F({\phi_{\pm}^{\textrm{coh}}})$ with respect to TMEPADFS (TMEPSDFS) parameter
$\alpha$ and for different values of $n$ and $k$ with $f(n)=1$.

In order to find the teleportation fidelity using the $f$-deformed two-mode
entangled resources, coefficients $C_{m}^{\pm}(n,k,\alpha)$ in (III.1) will be
replaced by $\mathbb{C}_{m}^{\pm}(n,k,\alpha)$. It can be seen that the
fidelity depends only on the Fock parameter $n$, photon-addition/subtraction
number $k$, coherent state amplitude $\alpha$ of the entangled resources and
is free from the amplitude $\alpha_{0}$ of the input coherent state.

The fidelity for teleporting a coherent state is plotted as a function of
amplitude $\alpha$ of entangled resources in Fig. 1. It is noted that for
$n=k=0$, the TMEPADFS switches to a two-mode entangled coherent state and the
corresponding fidelity approaches the maximum threshold limit of 0.5 (see Fig.
1(a)). The fidelity with the TMEPADFS is less than that with the original no
photon-added DFS. It is also observed that TMEPADFS yields lower fidelity
value compared to TMEPSDFS although PADFS is more nonclassical than PSDFS [1,
9]. Moreover, the fidelity plot goes up with an increase in the photon
subtraction number but diminishes as the photon addition number is increasing.
For a fixed value of photon-addition/subtraction number $k$, an increase in
the Fock parameter $n$ leads to weak fidelity, while varying the displacement
parameter $\alpha$ from 0 to 3 increases the fidelity value in the beginning
followed by a collapse. In case of $n=0$, PSDFS reduces to the elementary
coherent state, causing fidelity to coincide with that of DFS for different
values of the photon-subtraction number $k$ (see Fig. 1(d)). We see that for
the linear intensity-dependent function $f(n)=1$, the fidelity remains below
the classical limit of 1/2, thus fails to prove linear TMEPADFS and TMEPSDFS
relevant for quantum teleportation.

To overcome this challenge, we consider next the teleportation fidelity of a
coherent state with different $f$-deformed non-Gaussian entangled resources.
The fidelities corresponding to the nonlinear intensity-dependent functions
$f(n)=\sqrt{n}$ and $f(n)=1/\sqrt{n}$ are plotted in Figs. 2 and 3,
respectively. Comparing with panels I and II of Fig. 1, we can see that the
qualitative behaviours are very similar to the previous cases of linear
TMEPADFS and TMEPSDFS.

Figure 2: (Color online) Behaviour of fidelity
$F({\varphi_{\pm}^{\textrm{coh}}})$ associated to $f$-deformed entangled
resource states, plotted with respect to $\alpha$ and for different values of
$n$ and $k$ with $f(n)=\sqrt{n}$. Figure 3: (Color online) Plot of fidelity
$F({\varphi_{\pm}^{\textrm{coh}}})$ as a function of $\alpha$ for different
values of $n$, $k$ and $f(n)=1/\sqrt{n}$.

In Fig. 2, the fidelity initially experiences a decline and then ascends with
the non-Gaussian state channel parameter $\alpha$. However, it also remains
unable to cross the classical boundary. As the Fock parameter or photon
addition number increases, fidelity decreases while an increase in the photon
subtraction parameter results enhancement in fidelity. Fig. 3 shows a
different fidelity pattern for the intensity-dependent function
$f(n)=1/\sqrt{n}$. In this case, the fidelity plot outperforms the threshold
limit that leads to a significant advantage of using $f$-deformed resource
states. Similar to Fig. 2, an increase in the Fock parameter or photon
addition number corresponds to a fidelity reduction, while an increase in the
photon subtraction parameter leads to improvement of fidelity.

###  III.2 Input Squeezed State

Once the characteristic functions of the two-mode entangled photon-
added/subtracted displaced Fock states are derived, we now proceed to
calculate the fidelity for teleporting an input squeezed vacuum state
$\ket{z}$ with $z=r\,e^{i\zeta}$. The BK protocol is followed for teleporting
an input quantum state between two distant physical systems. The analytic
expression of fidelity using TMEPADFS (TMEPSDFS) channel turns out to be

$\displaystyle F({\phi_{\pm}^{\textrm{squ}}})$
$\displaystyle=\frac{1}{2\cosh{r}}\sum_{m,l=0}^{\infty}\frac{(C_{m}^{\pm}(n,k,\alpha))^{*}C_{l}(n,k,\alpha)}{2^{(m+l)/2}}\sqrt{\frac{l!}{m!}}$
$\displaystyle\times(-1)^{m-l}\sum_{j,k=0}^{m,l}{m\choose
j}\sum_{r,s=0}^{k,l-k}{j\choose k-r}{m-j\choose l-k-s}\frac{(-1)^{r+s}}{r!s!}$
$\displaystyle\times\sum_{t=0}^{j-k+r+s}\binom{j-k+r+s}{t}\frac{1}{2^{t}}\left(\frac{\tanh{r}}{2}\right)^{r+s-t+\frac{m-l}{2}}$
$\displaystyle\times\frac{(k-j+r+s+m-l)!}{(k-j+r+s+m-l-t)!!}$
$\displaystyle\times e^{i\zeta{(k-j+m-l)/2}}(j-k+r+s-t-1)!!$ (15)

if $k-j+r+s+m-l-t$ and $j-k+r+s-t$ are even and

$\displaystyle n!!=\left\\{\begin{array}[]{lll}&n(n-2)(n-4)\ldots 4.2&\mbox{if
$n$ is even},\\\ \\\ &n(n-2)(n-4)\ldots 3.1&\mbox{if $n$ is
odd}\end{array}\right.$

$F({\phi_{\pm}^{\textrm{squ}}})$ are plotted as a function of $\alpha$ in Fig.
4.

Figure 4: (Color online) Fidelity of teleportation
$F({\phi_{\pm}^{\textrm{squ}}})$ as a function of $\alpha$ with $r=1$ and
$f(n)=1$.

We notice that in all the cases under consideration,
$F({\phi_{+}^{\textrm{squ}}})$ initially increases until reaching a maximum
limit and subsequently decreases. The fidelity pattern for an input squeezed
state and TMEPADFS/TMEPSDFS/TMEDFS channels matches that of an input coherent
state, although a bit lower. Notably, it equals to the fidelity of the input
coherent state when the squeezing parameter $r$ is zero. To explore the
fidelity for nonlinear intensity-dependent functions, we have plotted this for
teleporting an input squeezed state via $f$-deformed channels such as
TMEPADFS/TMEPSDFS/TMEDFS with $f(n)=\sqrt{n}$ and $f(n)=1/\sqrt{n}$ in Figs. 5
and 6, respectively.

Figure 5: (Color online) Fidelity $F({\varphi_{\pm}^{\textrm{squ}}})$ is
plotted with respect to $\alpha$ with $r=1$ and $f(n)=\sqrt{n}$. Figure 6:
(Color online) Variation of fidelity $F({\varphi_{\pm}^{\textrm{squ}}})$ with
respect to $\alpha$ with $r=1$ and $f(n)=1/\sqrt{n}$.

Fig. 5 clearly shows that for the intensity-dependent function
$f(n)=\sqrt{n}$, the teleportation fidelity initially decreases and then
increases with the non-Gaussian state channel parameter $\alpha$. However, it
remains unable to surpass the threshold limit. Increasing the Fock parameter
or photon addition number results in fidelity reduction, while an increase in
the photon subtraction number leads to fidelity improvement. Similarly in Fig.
6, we observe that for the intensity-dependent function $f(n)=1/\sqrt{n}$, the
teleportation fidelity first increases and then decreases with the non-
Gaussian channel parameter $\alpha$. Like the previous case, it fails to cross
the threshold limit. The fidelity decreases with an increase in the Fock
parameter or photon addition number, while increases with an increase in the
photon subtraction number.

To examine the performance of TMEPADFS/TMEPSDFS channels for teleporting an
input squeezed vacuum state with respect to the squeezing parameter $z$, the
variations of $F({\phi_{\pm}^{\textrm{squ}}})$ and
$F({\varphi_{\pm}^{\textrm{squ}}})$ are plotted with respect to $z$ for
$f(n)=1$, $\sqrt{n}$ and $1/\sqrt{n}$ in Figs. 7, 8 and 9, respectively.

Figure 7: (Color online) Fidelity $F({\phi_{\pm}^{\textrm{squ}}})$ as a
function of squeezing parameter $z$ with $\alpha=1$ and $f(n)=1$. Figure 8:
(Color online) Fidelity $F({\varphi_{\pm}^{\textrm{squ}}})$ with respect to
squeezing parameter $z$ of an input squeezed state $\ket{z}$ with $\alpha=1$
and $f(n)=\sqrt{n}$. Figure 9: (Color online) Variation of fidelity
$F({\varphi_{\pm}^{\textrm{squ}}})$ as a function of $z$ with $\alpha=1$ and
$f(n)=1/\sqrt{n}$.

These figures exhibit that the fidelity remains below the classical limit for
different values of $f(n)$. Thus for squeezed state inputs, using nonlinear
intensity-dependent functions has not been proved beneficial in terms of the
fidelity. It is also evident that the fidelity reaches its maximum at a low
squeezing regime and subsequently decreases.

##  IV Entanglement and EPR Correlation of $f$-deformed States

Entangled resources are highly advantageous in various quantum information
processing tasks such as continuous-variable quantum teleportation [14],
quantum dense coding [31], and entanglement swapping [32]. Entanglement is a
quantum mechanical phenomenon where two or more particles become correlated in
such a way that the state of one particle cannot be described independently of
the state of the other particle(s) [33]. In other words, the particles become
“entangled” and their properties are reliant with each other. For any
bipartite pure state, the entanglement can be quantified by using the partial
von Neumann entropy as $E(\Psi_{12})=S(\rho_{1})=-\Tr[\rho_{1}\ln(\rho_{1})]$
i.e. entanglement of formation [34]. The entanglement for $f$-deformed
TMEPADFS and TMEPSDFS can be calculated as

$\displaystyle E(\varphi_{\pm})$
$\displaystyle=-\sum_{m,l,j=0}^{\infty,\infty,m}\frac{|\mathbb{C}_{m}^{\pm}(n,\alpha)\mathbb{C}_{l}^{\pm}(n,\alpha)|}{2^{(m+l)/2}}\sqrt{{m\choose
j}{l\choose l-m+j}}$
$\displaystyle\times\ln(\frac{|\mathbb{C}_{m}^{\pm}(n,\alpha)\mathbb{C}_{l}^{\pm}(n,\alpha)|}{2^{(m+l)/2}}\sqrt{{m\choose
j}{l\choose l-m+j}})$ (17)

The entanglement of these states is plotted in Figs. 10, 11 and 12.

Figure 10: (Color online) Entanglement of formation $E(\phi_{\pm})$ as a
function of TMEPADFS/TMEPSDFS parameter $\alpha$ for different values of $n$
and $k$, $f(n)=1$. Figure 11: (Color online) Entanglement of formation
$E(\varphi_{\pm})$ as a function of TMEPADFS/TMEPSDFS parameter $\alpha$ for
different values of $n$ and $k$ and for $f(n)=\sqrt{n}$. Figure 12: (Color
online) Entanglement of formation $E(\varphi_{\pm})$ as a function of
TMEPADFS/TMEPSDFS parameter $\alpha$ for different values of $n$ and $k$ and
for $f(n)=1/\sqrt{n}$.

We observe that the entanglement of all the states under consideration first
increases upto the maximum value and then decreases rapidly with respect to
the displacement parameter $\alpha$. Further, increasing the photon addition
enhances the entanglement while increasing the number of subtracted photons
reduces the entanglement value. Also adding photons seems to be more
significant than subtracting. Increasing the Fock state parameter $n$
increases the amount of entanglement.

In addition to the degree of entanglement, entangled states can also be
characterized by the Einstein-Podolsky-Rosen (EPR) correlations between phase-
quadrature components of the two modes. Analogous to the position and momentum
operators of a macroscopic particle, the phase-quadrature operators of each
mode are defined as $x_{j}=\frac{1}{\sqrt{2}}(a_{j}+a_{j}^{\dagger})$ and
$p_{j}=\frac{1}{i\sqrt{2}}(a_{j}-a_{j}^{\dagger})$ for $j=1,2$. In the vacuum
state, both the variances $\Delta(x_{1}-x_{2})^{2}$ and
$\Delta(p_{1}+p_{2})^{2}$ are equal to one. For any two-mode classical state,
these variances $\Delta(x_{1}-x_{2})^{2}$ and $\Delta(p_{1}+p_{2})^{2}$ are
greater than 1. In the EPR state,
$\Delta(x_{1}-x_{2})^{2}=\Delta(p_{1}+p_{2})^{2}=0$. That means $x_{1}$ and
$p_{1}$ of the first mode can be exactly estimated by measured results of
$x_{2}$ and $p_{2}$ of the second mode or vice versa. In this sense, the ideal
EPR correlation exists between the two modes. There may be some two-mode
states, e.g. two-mode squeezed state that possess the EPR correlation beyond
the limit of the vacuum state, i.e., both the variances
$\Delta(x_{1}-x_{2})^{2}$ and $\Delta(p_{1}+p_{2})^{2}$ are less than 1. This
quantum correlation is a key ingredient for realizing quantum teleportation of
continuous variables.

For all the states, $\Delta(x_{1}-x_{2})^{2}=\Delta(p_{1}+p_{2})^{2}$. In
Figs. 13, 14 and 15, the variance $\Delta(x_{1}-x_{2})^{2}$ for TMEDFS,
TMEPADFS, TMEPSDFS is plotted.

Figure 13: (Color online) Variance $\Delta(x_{1}-x_{2})^{2}$ for the state
$\ket{\phi_{\pm}}$ with $f(n)=1$ and as a function of $\alpha$ for different
values of $n$ and $k$. Figure 14: (Color online) Variance
$\Delta(x_{1}-x_{2})^{2}$ for the state $\ket{\varphi_{\pm}}$ with
$f(n)=\sqrt{n}$ and as a function of $\alpha$ for different values of $n$ and
$k$. Figure 15: (Color online) Variance $\Delta(x_{1}-x_{2})^{2}$ for the
state $\ket{\varphi_{\pm}}$ with $f(n)=1/\sqrt{n}$ and as a function of
$\alpha$ for different values of $n$ and $k$.

It has been found that the variance $\Delta(x_{1}-x_{2})^{2}$ decreases
rapidly with respect to $\alpha$. The EPR correlation for TMEPADFS as well as
TMEPSDFS decreases as $n$ changes from 0 to 2. We have further noticed that
the variance of two-mode entangled photon-added displaced Fock state is less
than that of two-mode entangled photon-subtracted displaced Fock state.
Therefore, one may expect that the quality of quantum teleportation of
continuous-variables can be improved by using two-mode entangled photon-
subtracted displaced Fock states.

##  V Conclusion

In this work, the success probability for employing $f$-deformed two-mode
entangled resources in continuous-variable quantum teleportation is evaluated.
The entanglement properties of the non-Gaussian states are also investigated
which are constructed by mixing a single-mode PADFS/PSDFS and a vacuum state
at a 50:50 beam-splitter. We study quantum teleportation of Braunstein and
Kimble protocol for coherent states and squeezed states, with the resulting
non-Gaussian states as entangled resources. For coherent states, the
analytical expression of fidelity is derived. It is noted that the fidelity is
independent of the amplitude of coherent states to be teleported. For the
squeezed vacuum input, the fidelity is state dependent. For all the states to
be teleported, we have observed that the fidelity with the TMEPSDFS is higher
than that with the TMEPADFS and it decreases as Fock parameter $n$ increases.
The fidelity for teleporting coherent state with the non-Gaussian $f$-deformed
two-mode entangled resources moves across the benchmark value of 1/2.
Moreover, in the weak squeezing regime of the original squeezed state
$\ket{z}$, the TMEPADFS and TMEPSDFS can lead to the highest fidelity,
although fail to cross the classical limit.

We have also observed that in a particular region, the partial von Neumann
entropy of all the resulting states is greater than that of the original two-
mode entangled coherent state in a particular range of $\alpha$. The photon
addition appears to be a more efficient operation for enhancing the amount of
entanglement. As the EPR correlation between phase-quadrature components of
the two modes in the bipartite states signals the existence of entanglement,
we have found that TMEPSDFS states have stronger EPR correlation than the
TMEPADFS. Moreover, in the low-energy regime of the input squeezed vacuum
state, the plot for variance collapses quickly.

In addition, if the nonlinear intensity-dependent function can be fixed as
$f(n)=1/n^{p}$, then for a large positive real number $p$, it becomes possible
to achieve teleportation fidelity close to 1 for any of the input states.

ACKNOWLEDGEMENT

Deepak’s work is supported by the Council of Scientific and Industrial
Research (CSIR), Govt. of India (Award no. 09/1256(0006)/2019-EMR-1). A. C.
acknowledges DST SERB for the support provided through the project number
SUR/2022/000899.

## References

  * Malpani  _et al._ [2023] P. Malpani, K. Thapliyal, J. Banerji, and A. Pathak, Enhancement of non-Gaussianity and nonclassicality of photon-added displaced Fock state: A quantitative approach, Annalen der Physik 535, 2200261 (2023).
  * Deepak and Chatterjee [2023] Deepak and A. Chatterjee, Nonclassicality versus quantum non-Gaussianity of photon-subtracted displaced Fock state, Canadian Journal of Physics 101, 560 (2023).
  * Bartlett  _et al._ [2002] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Efficient classical simulation of continuous variable quantum information processes, in  _Quantum Information with Continuous Variables_ (Springer Netherlands, 2002) pp. 47–55.
  * Bartlett and Sanders [2002] S. D. Bartlett and B. C. Sanders, Efficient classical simulation of optical quantum information circuits, Physical Review Letters 89, 207903 (2002).
  * Eisert  _et al._ [2002] J. Eisert, S. Scheel, and M. B. Plenio, Distilling Gaussian states with Gaussian operations is impossible, Physical Review Letters 89, 137903 (2002).
  * Fiurášek [2002] J. Fiurášek, Gaussian transformations and distillation of entangled Gaussian states, Physical Review Letters 89, 137904 (2002).
  * Giedke and Ignacio Cirac [2002] G. Giedke and J. Ignacio Cirac, Characterization of Gaussian operations and distillation of Gaussian states, Physical Review A 66, 032316 (2002).
  * Niset  _et al._ [2009] J. Niset, J. Fiurášek, and N. J. Cerf, No-go theorem for Gaussian quantum error correction, Physical Review Letters 102, 120501 (2009).
  * Malpani  _et al._ [2019] P. Malpani, N. Alam, K. Thapliyal, A. Pathak, V. Narayanan, and S. Banerjee, Lower‐ and higher‐order nonclassical properties of photon added and subtracted displaced Fock states, Annalen der Physik 531, 10.1002/andp.201800318 (2019).
  * Bennett  _et al._ [1993] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels, Physical Review Letters 70, 1895 (1993).
  * Bouwmeester  _et al._ [1997] D. Bouwmeester, J.-W. Pan, K. Mattle, M. Eibl, H. Weinfurter, and A. Zeilinger, Experimental quantum teleportation, Nature 390, 575 (1997).
  * Boschi  _et al._ [1998] D. Boschi, S. Branca, F. De Martini, L. Hardy, and S. Popescu, Experimental realization of teleporting an unknown pure quantum state via dual classical and Einstein-Podolsky-Rosen channels, Physical Review Letters 80, 1121 (1998).
  * Vaidman [1994] L. Vaidman, Teleportation of quantum states, Physical Review A 49, 1473 (1994).
  * Braunstein and Kimble [1998] S. L. Braunstein and H. J. Kimble, Teleportation of continuous quantum variables, Physical Review Letters 80, 869 (1998).
  * Bowen  _et al._ [2003] W. P. Bowen, N. Treps, B. C. Buchler, R. Schnabel, T. C. Ralph, H.-A. Bachor, T. Symul, and P. K. Lam, Experimental investigation of continuous-variable quantum teleportation, Physical Review A 67, 032302 (2003).
  * Furusawa  _et al._ [1998] A. Furusawa, J. L. Sørensen, S. L. Braunstein, C. A. Fuchs, H. J. Kimble, and E. S. Polzik, Unconditional quantum teleportation, Science 282, 706 (1998).
  * Takei  _et al._ [2005] N. Takei, H. Yonezawa, T. Aoki, and A. Furusawa, High-fidelity teleportation beyond the no-cloning limit and entanglement swapping for continuous variables, Physical Review Letters 94, 220502 (2005).
  * van Enk [1999] S. J. van Enk, Discrete formulation of teleportation of continuous variables, Physical Review A 60, 5095 (1999).
  * Vukics  _et al._ [2002] A. Vukics, J. Janszky, and T. Kobayashi, Nonideal teleportation in coherent-state basis, Physical Review A 66, 023809 (2002).
  * Hofmann  _et al._ [2000] H. Hofmann, T. Ide, T. Kobayashi, and A. Furusawa, Fidelity and information in the quantum teleportation of continuous variables, Physical Review A 62, 062304 (2000).
  * Navarrete-Benlloch [2015] C. Navarrete-Benlloch, Quantum information with continuous variables, in  _An Introduction to the Formalism of Quantum Information with Continuous Variables_ (Morgan & Claypool Publishers, 2015) pp. 4–1–4–57.
  * Deepak and Chatterjee [2022] Deepak and A. Chatterjee, Realistic continuous-variable quantum teleportation using a displaced Fock state channel, Quantum Information Processing 21, 10.1007/s11128-022-03484-y (2022).
  * Rosen [1979] N. Rosen, Can quantum-mechanical description of physical reality be considered complete?, in  _Albert Einstein_ (Vieweg+Teubner Verlag, 1979) pp. 57–67.
  * Man’ko  _et al._ [1997] V. I. Man’ko, G. Marmo, E. C. G. Sudarshan, and F. Zaccaria, $f-$oscillators and nonlinear coherent states, Physica Scripta 55, 528 (1997).
  * Faghihi [2020] M. J. Faghihi, Generalized photon added and subtracted $f-$deformed displaced Fock states, Annalen der Physik 532, 2000215 (2020).
  * Marian and Marian [2006] P. Marian and T. A. Marian, Continuous-variable teleportation in the characteristic-function description, Physical Review A 74, 042306 (2006).
  * Quintero and Ladera [2011] W. Quintero and C. L. Ladera, Teleportation of displaced Fock states: Fidelity and their teleported photon number distributions, Journal of Physics: Conference Series 274, 012146 (2011).
  * Dell’Anno  _et al._ [2007] F. Dell’Anno, S. De Siena, L. Albano, and F. Illuminati, Continuous-variable quantum teleportation with non-Gaussian resources, Physical Review A 76, 022301 (2007).
  * Dell’Anno  _et al._ [2010] F. Dell’Anno, S. De Siena, and F. Illuminati, Realistic continuous-variable quantum teleportation with non-Gaussian resources, Physical Review A 81, 012333 (2010).
  * Kumar and Arora [2023] C. Kumar and S. Arora, Success probability and performance optimization in non-Gaussian continuous-variable quantum teleportation, Phys. Rev. A 107, 012418 (2023).
  * Braunstein and Kimble [2000] S. L. Braunstein and H. J. Kimble, Dense coding for continuous variables, Phys. Rev. A 61, 042302 (2000).
  * Tan [1999] S. M. Tan, Confirming entanglement in continuous variable quantum teleportation, Physical Review A 60, 2752 (1999).
  * Yang and Li [2009] Y. Yang and F.-L. Li, Entanglement properties of non-Gaussian resources generated via photon subtraction and addition and continuous-variable quantum-teleportation improvement, Physical Review A 80, 022315 (2009).
  * Bennett  _et al._ [1996] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Mixed-state entanglement and quantum error correction, Physical Review A 54, 3824 (1996).


# Modeling ice crystal growth using the lattice Boltzmann method

Q. Tan  S.A. Hosseini<EMAIL_ADDRESS>A. Seidel-Morgenstern  D.
Thévenin  H. Lorenz  Laboratory of Fluid Dynamics and Technical Flows,
University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany
Department of Mechanical and Process Engineering, ETH Zürich, 8092 Zürich,
Switzerland  Max Planck Institute for Dynamics of Complex Technical Systems
(MPI DCTS), 39106 Magdeburg, Germany

###### Abstract

Given the multitude of growth habits, pronounced sensitivity to ambient
conditions and wide range of scales involved, snowflake crystals are one of
the most challenging systems to model. The present work focuses on the
development and validation of a coupled flow/species/phase solver based on the
lattice Boltzmann method. It is first shown that the model is able to
correctly capture species and phase growth coupling. Furthermore, through a
study of crystal growth subject to ventilation effects, it is shown that the
model correctly captures hydrodynamics-induced asymmetrical growth. The
validated solver is then used to model snowflake growth under different
ambient conditions with respect to humidity and temperature in the plate-
growth regime section of the Nakaya diagram. The resulting crystal habits are
compared to both numerical and experimental reference data available in the
literature. The overall agreement with experimental data shows that the
proposed algorithm correctly captures both the crystal shape and the onset of
primary and secondary branching instabilities. As a final part of the study
the effects of forced convection on snowflake growth are studied. It is shown,
in agreement with observations in the literature, that under such condition
the crystal exhibits non-symmetrical growth. The non-uniform humidity around
the crystal due to forced convection can even result in the coexistence of
different growth modes on different sides of the same crystal.

###### keywords:

phase-field models, Lattice Boltzmann method, anisotropic snowflake growth,
Ventilation effects.

###### MSC:

[2010] 00-01, 99-00

††journal: …

##  1 Introduction

Ice crystals and their growth are of interest in many fields such as
environmental sciences, agriculture and industries as aviation (encountering
de-icing issues). Due to the wide spectrum of habits they take on, they have
been the topic of scientific research for decades [1, 2, 3]. Johannes Kepler
was the first person to explore the growth mechanisms of snowflakes in 1611,
and attempted to explain the possible origins of snow crystal symmetry [4].
With the development of photography in the late 19th century, Wilson Bentley
[5] collected in 1931 thousands of snow crystal images [6]. In 1938, Ukichiro
Nakaya and his co-workers [7, 8, 9] conducted comprehensive experimental
studies of ice crystal growth to determine the relationship between growth
conditions and crystal shape. The different growth modes and associated habits
as a function of temperature and supersaturation (humidity) are illustrated in
Fig. 1.

Figure 1: The snow crystal morphology diagram, showing the morphology of ice
crystals growing from water vapor in air at 1 bar as a function of temperature
and supersaturation. This figure is taken from [10, 11].

Overall, based on the dominant growth direction, snow-flakes can be classified
as pertaining to one of these two categories: plates, or columns/needles. The
growth mode is essentially dictated by the temperature, while growth rate, and
the associated instabilities, are affected by both temperature and humidity.
As observed in Fig. 1, crystal habits alternate between flat and columnar.
Transitions between the different modes, i.e. plate and column, occur at
around -5, -10 and -20∘C. Furthermore, higher humidity and thus, higher
supersaturation as crystallization driving force contributes to more
pronounced instability effects and therefore more complex shapes. Despite the
efforts by Ukichiro Nakaya [9] and other researchers, much of the
phenomenology behind the growth of snowflakes still remains uncertain [12, 13,
14, 15]. In 1958, Mason [16] suggested that the basic habit is determined by
the surface diffusion rates. Later, Nelson and Knight [17] suggested that
layer nucleation rates influence the morphology of snow crystals. Libbrecht
[14] suggested that the different crystal habits and instabilities appearing
during growth are direct consequences of two main competing mechanisms: vapor
diffusion in air and kinetics of the water molecules attachment to the crystal
surface [8, 18, 19]. The geometrically ordered potential field on the crystal
surface dictates the growth into faceted structures (resulting from the
molecules’ arrangement within the crystal lattice), while diffusion
contributes to growth instabilities that produce dendritic branching and
associated complex structures. Attempts to further explain the growth
mechanism of snowflakes are still ongoing.

While most of the previous research was focused on experimental studies of ice
growth, widespread effort has also been put on developing mathematical models
and numerical tools. Microscopic models and simulations such as Molecular
Dynamic simulations [20] have been conducted to analyze fundamental processes
such as surface adsorption and diffusion. While physically sound, these
simulations are limited in time and space and are not efficient for
simulations of a full snowflake (especially at later stages when instabilities
appear and the crystal grow in size). So-called mesoscopic formulations are an
interesting alternative. The cellular automaton of Gravner and Griffeath [21]
is an illustration of such approaches. While providing spectacular results in
predicting faceted growth, such models lack established connections with
physical processes and parameters [22]. Furthermore, reliable models
introducing additional physics such as interaction with a flow field are yet
to be developed [23]. At the macroscopic level, crystal growth is modeled
using either sharp or diffuse interface formulations. Various snowflake
morphologies were simulated by Barrett et al. [24] with a sharp interface
model. However, only small supersaturations could be considered because of the
numerical cost of the interface parametrization [25]. In recent years, the
phase-field model [26], a type of diffuse interface formulation, has become
one of the most popular methods for the simulation of crystal growth. The
phase-field model is a powerful tool to simulate interface development in the
crystallization process as the model does not require explicit tracking of the
interface via front-tracking algorithms. In addition, the non-linear partial
differential equations are obtained from the principles of non-equilibrium
thermodynamics making the interface dynamics consistent without the need for
explicit boundary treatments. Over the past several decades, few successful
attempts are reported to model faceted snowflake growth [27, 28]. Comparisons
with experimental data were rather promising. lbmLBMlattice Boltzmann method
lbLBlattice Boltzmann nsNSNavier-Stokes pdePDEpartial differential equation

The , developed over the past decades, has become a popular alternative to
classical solvers for the equations [29]. It has since been extended to a
variety of applications and flow regimes such as multi-phase flows [30], flows
in porous media, turbulence and multi-component flows [31, 32, 33]. It has
also been used, in combination with classical solvers for the solid phase, to
simulate crystal growth [23, 34, 35, 36, 37, 38, 39, 40]. While initially
developed as a discrete solver for the Boltzmann equation in the hydrodynamic
regime, it has also widely been used as a solver for different types of
parabolic via appropriate parametrization of the equilibrium state and the
collision-streaming operators. In that same spirit a number of work have
proposed -based formulations to solve the phase-field evolution equation.
Younsi et al. [41] modified the standard equations and equilibrium
distribution function to introduce anisotropic surface tension and growth rate
effects. The proposed model was successfully used to simulate anisotropic and
dendritic growth. While readily applied to generic systems, these models have
rarely been used to simulate realistic systems such as snowflake growth.

The aim of the present work is to present a pure model able to correctly
capture not only different growth modes of snowflakes in the platelet regime,
but also effects caused by the surrounding convective field. For this purpose,
after an overview of the approach proposed to model crystal growth, it will
first be validated against generic test-cases well-documented in the
literature (both with and without ventilation effects). The solver is then
used to model the growth of a single ice crystal under different conditions
regarding temperature and supersaturation. The obtained shapes are validated
(qualitatively) against both experimental and numerical data available in the
literature. It will also be shown that the growth modes are in agreement with
those predicted in the Nakaya diagram. Finally, the effect of fluid flow on
dendritic growth is studied. Interestingly, rarely studied growth habits are
observed in the simulations;For example triangular shapes as observed only
when crystals are exposed to a background convection.

##  2 Theoretical background

While influenced by both temperature and humidity (referred to as
supersaturation in what follows), given the more pronounced effect of the
latter on faceted growth and dendritic instabilities, the temperature field
will be considered to be uniform during the present study. As such, the phase-
field growth dynamics will involve two coupled equations (apart from the flow
field solver). One will describe the morphology of the snowflakes through a
phase parameter $\phi$; the second one is an advection/diffusion-type equation
for the supersaturation field $U$.

###  2.1 Diffuse-interface formulation: governing equations

The kinetics of phase-field are described by a parameter $\phi$ designating
solid and fluid phases respectively when reaching +1 (ice) and -1 (vapor) and
the reduced supersaturation of water vapor defined by
$U=(c-c_{sat}^{S})/c_{sat}^{S}$, where $c_{sat}^{S}(T)$ is the saturation
number density of vapor of ice at temperature $T$. The space/time evolution
equations are written as [27, 41, 42, 43]:

$\tau_{0}a_{s}^{2}(\textbf{n})\frac{\partial\phi}{\partial t}=\\\
W_{0}^{2}\bm{\nabla}\cdot\left(a_{s}^{2}(\textbf{n})\right)\bm{\nabla}\phi+W_{0}^{2}\bm{\nabla}\cdot\left(|\bm{\nabla}\phi|^{2}\frac{\partial[a(\textbf{n})^{2}]}{\partial\bm{\nabla}\phi}\right)\\\
+(\phi-\phi^{3})+\lambda U(1-\phi^{2})^{2},$ (1)

and:

$\frac{\partial U}{\partial
t}+\left(\frac{1-\phi}{2}\right)\bm{u}\cdot\bm{\nabla}U=D\bm{\nabla}\cdot\left(q(\phi)\bm{\nabla}U\right)-\frac{L_{sat}}{2}\frac{\partial\phi}{\partial
t},$ (2)

defining for the later analysis $\tau=\tau_{0}a_{s}^{2}(\textbf{n})$. Here,
$\tau_{0}$ denotes the characteristic time and $W_{0}$ denotes the
characteristic width of the diffuse interfaces. In Eq. 1, the coefficient
$\lambda=\frac{15L^{2}}{16Hc_{p}T_{m}}$ describes the strength of the coupling
between the phase-field and the species field. $1/\lambda$ is a dimensionless
measure of the barrier height $H$ of the double-well potential. The latent
heat of melting is shown as $L$. The specific heat capacity $c_{p}$ is assumed
to be the same in the two phases (symmetric model).
$\textbf{n}=-\frac{\bm{\nabla}\phi}{\left|\bm{\nabla}\phi\right|}$ is the unit
vector normal to the crystal interface –pointing from solid to the fluid.
$a_{s}(\textbf{n})$ is the surface tension anisotropy function which here, in
the context of the snowflake growth studies, is defined as:

$a_{s}(\textbf{n})=1+\epsilon_{xy}\cos(6\theta),$ (3)

rhsRHSright hand side lhsLHSleft hand side where
$\theta=\arctan(n_{y}/n_{x})$. $\epsilon_{xy}$ is a numerical parameter
characterizing the anisotropy strength. The second and third terms on the of
Eq. 1 control the anisotropic growth of dendrites in snowflakes.
$\phi-\phi^{3}$ is the derivative of the double-well potential. The last term
in Eq. 1 is a source term accounting for the coupling between the
supersaturation $U$ and the order parameter $\phi$. $(1-\phi^{2})^{2}$ is an
interpolation function minimizing the bulk potential at $\phi=\pm 1$.
In Eq. 2, $\bm{u}$ denotes the fluid local velocity while $q(\phi)=(1-\phi)$
is a function canceling out diffusion within the solid. As such solute
transport is assumed to take place only within the fluid phase (one-side
model). The parameter $D$ is the diffusion coefficient of water vapor in air
as applied in Section 3.2.1. $L_{sat}$ describes the absorption rate of water
vapor at the interface of the snow crystals.
The thermal capillary length is defined as:

$d_{0}=\frac{\gamma_{0}T_{M}c_{p}}{L^{2}},$ (4)

where the excess isotropic free energy of the solid-liquid interface
$\gamma_{0}$ is given by:

$\gamma_{0}=IW_{0}H,$ (5)

where $W_{0}$ is the interface thickness, $H$ is the barrier height of the
double-well potential and $I=2\sqrt{2}/3$. Hence, the matched asymptotic
expansions provide a relation between the phase-field and sharp-interface
parameters given by:

$d_{0}=a_{1}\frac{W_{0}}{\lambda}.$ (6)

The expression for the kinetic coefficient is:

$\beta=a_{1}\left(\frac{\tau_{0}}{W_{0}\lambda}-a_{2}\frac{W_{0}}{D}\right)$
(7)

where $a_{2}$ is a constant and depends on the choice of the function
$g(\phi)$. In this paper, the standard form
$g(\phi)=\phi-2\phi^{3}/3+\phi^{5}/5$ will be used. For the model,
$a_{1}\approx 0.8839$ and $a_{2}\approx 0.6267$. These relations make it
possible to choose phase-field parameters for prescribed values of the
capillary length (surface energy) and the interface mobility (interface
kinetic coefficient). Note that the interface width $W_{0}$ is a parameter
that can be freely chosen in this formulation; the asymptotic analysis remains
valid as long as $W_{0}$ remains much smaller than any length scale present in
the sharp-interface solution of the considered problem [26, 44, 45].
It follows from Eq.(7) that $\beta$ can be made to vanish provided that $\tau$
is chosen to be:

$\tau=\tau_{0}\left[A(n)^{2}\right].$ (8)

Therefore,

$\tau_{0}=a_{2}\lambda\frac{W_{0}^{2}}{D}.$ (9)

The results should be independent of $\lambda$ when they are converged.
Decreasing $\lambda$ corresponds physically to decreasing the interface width
while increasing at the same time the height of the double-well potential so
as to keep the surface energy and hence $d_{0}$ fixed.
The physical dimensions of $W_{0}$, $\tau_{0}$ and $\lambda$ are respectively
$[W_{0}]\equiv[\mathcal{L}]$, $[\tau_{0}]\equiv[\mathcal{T}]$,
$[\lambda]\equiv[-]$, $[D]\equiv[\mathcal{L}^{2}/\mathcal{T}]$ where
$[\mathcal{L}]$ indicates the length dimension and $[\mathcal{T}]$ indicates
the time dimension.

###  2.2 Lattice Boltzmann formulation

####  2.2.1 Flow field solver

The target system of equations describing the flow field behavior (i.e.
incompressible -continuity equations) are modeled using the classical
isothermal formulation consisting of the now-famous stream-collide operators:

$f_{\alpha}\left(\bm{x}+\bm{c}_{\alpha}\delta_{t},t+\delta_{t}\right)-f_{\alpha}\left(\bm{x},t\right)=\delta_{t}\Omega_{\alpha}\left(\bm{x},t\right)+\delta_{t}\bm{F},$
(10)

where $\bm{F}$ is the external force. Here, the external force $\bm{F}$
modeling interaction with the solid phase is given as [43]:

$F=-\frac{h\eta_{f}(1+\phi)^{2}(1-\phi)\bm{u}}{4W_{0}^{2}}$ (11)

where $h$ is a dimensionless constant, here $h=2.757$. Due to the absence of
fluid velocity in snow crystals, the fluid velocity $\bm{u}$ is updated as:

$\bm{u^{*}}=\frac{(1-\phi)}{2}\bm{u},$ (12)

where the updated fluid velocity $\bm{u^{*}}$ is taken into the momentum
equation. The above friction term acts as a distributed momentum sink that
gradually forces the liquid velocity to zero as $\phi\rightarrow 1$. The
collision operator $\Omega_{\alpha}$ follows the linear Batnagar-Gross-Krook
approximation:

$\Omega_{\alpha}=\frac{1}{\tau}\left[f^{(eq)}_{\alpha}-f_{\alpha}\right],$
(13)

edfEDFequilibrium distribution function where $f_{\alpha}^{(eq)}$ is the
discrete isothermal defined as:

$f_{\alpha}^{(eq)}=\rho\sum_{i}\frac{1}{2c_{s}^{2}}a^{(eq)}_{i}(\bm{u}):\mathcal{H}_{i}(\bm{c}_{\alpha}),$
(14)

where $a^{(eq)}_{i}$ and $\mathcal{H}_{i}(\bm{c}_{\alpha})$ are the
corresponding multivariate Hermite coefficients and polynomials of order $i$,
with $c_{s}$ the lattice sound speed corresponding to the speed of sound at
the stencil reference temperature [46]. Further information on the expansion
along with detailed expressions of the can be found in [47, 48, 49]. The
relaxation time $\tau$ is tied to the fluid kinematic viscosity as:

$\tau=\frac{\nu}{c_{s}^{2}}+\frac{\delta_{t}}{2}.$ (15)

It must be noted that conserved variables, i.e. density and momentum are
defined as moments of the discrete distribution function:

$\rho=\sum_{\alpha}f_{\alpha},$ (16)
$\rho\bm{u}=\sum_{\alpha}\bm{c}_{\alpha}f_{\alpha}.$ (17)

####  2.2.2 Advection-diffusion-reaction solver for supersaturation field

The space/time evolution equation of the water vapour supersaturation field is
modeled using an advection-diffusion-reaction -based discrete kinetic
formulation defined as[50]:

$g_{\alpha}\left(\bm{x}+\bm{c}_{\alpha}\delta_{t},t+\delta_{t}\right)-g_{\alpha}\left(\bm{x},t\right)=\delta_{t}\Omega_{\alpha}\left(\bm{x},t\right)+\delta_{t}\dot{\omega}_{\alpha},$
(18)

where $\dot{\omega}_{\alpha}$ is the source term defined as:

$\dot{\omega}_{\alpha}=-w_{\alpha}\frac{L_{sat}}{2}\frac{\partial\phi}{\partial
t}.$ (19)

The collision operator $\Omega_{\alpha}$ for the supersaturation field is:

$\Omega_{\alpha}=\frac{1}{\tau_{g}}\left[g^{(eq)}_{\alpha}-g_{\alpha}\right].$
(20)

where $g_{\alpha}^{(eq)}$ is the defined as:

$g^{(eq)}_{\alpha}=\omega_{\alpha}U\left[1+\frac{\bm{c}_{\alpha}\cdot\bm{u}}{c_{s}^{2}}\right].$
(21)

The supersaturation is computed locally as the zeroth-order moment of
$g_{\alpha}$:

$U=\sum_{\alpha}g_{\alpha},$ (22)

and the relaxation coefficient is tied to the water vapour diffusion
coefficient [51] in air:

$\tau=\frac{Dq(\phi)}{c_{s}^{2}}+\frac{\delta_{t}}{2}.$ (23)

####  2.2.3 Solver for phase-field equation

The phase-field equation is modeled using a modified scheme defined as [52,
53]:

$a_{s}^{2}(\bm{n})h_{\alpha}(\bm{x}+\bm{c}_{\alpha}\delta x,t+\delta t)=\\\
h_{\alpha}(\bm{x},t)-\left(1-a_{s}^{2}(\bm{n})\right)h_{\alpha}(\bm{x}+\bm{c}_{\alpha}\delta
x,t)-\\\
\frac{1}{\eta_{\phi}(\bm{x},t)}\left[h_{\alpha}(\bm{x},t)-h_{\alpha}^{eq}(\bm{x},t)\right]+\omega_{\alpha}Q_{\phi}(\bm{x},t)\frac{\delta
t}{\tau_{0}},$ (24)

where the scalar function $Q_{\phi}$ is the source term of the phase-field
defined by:

$Q_{\alpha}=(\phi-\phi^{3})+\lambda U(1-\phi^{2})^{2},$ (25)

while the $h_{i}^{eq}$ is defined as:

$h_{\alpha}^{eq}=\omega_{\alpha}\left(\phi-\frac{1}{c_{s}^{2}}\bm{c}_{\alpha}\cdot\frac{W_{0}^{2}}{\tau_{0}}|\bm{\nabla}\phi|^{2}\frac{\partial(a_{s}(\bm{n})^{2}}{\partial\bm{\nabla}\phi}\frac{\delta
t}{\delta x}\right).$ (26)

The local value of the order parameter $\phi$ is calculated as:

$\phi=\sum_{\alpha}h_{\alpha},$ (27)

The relaxation time $\eta_{\phi}$ is a function of position and time and must
be updated at each time-step as:

$\eta_{\phi}=\frac{1}{c_{s}^{2}}a_{s}^{2}(\bm{n})\frac{W_{0}^{2}}{\tau_{0}}+\frac{\delta_{t}}{2},$
(28)

##  3 Simulations and numerical studies

For the sake of clarity, this part is organized into two different
subsections. In the first subsection results from simulations of generic
systems, intended as validation are presented and discussed. The results for
snowflakes in the plate regime are presented in the second subsection.

###  3.1 Validation studies

####  3.1.1 Effect of directional derivatives of gradients in LB scheme

In order to calculate the derivatives of $a_{s}(\bm{n})$ with respect to
$\partial_{x}\phi$ and $\partial_{y}\phi$ involved in the second term in Eq.
1, different choices are available. Overall, one can either evaluate the
gradient using classical finite-difference approximations, e.g. the central
difference second-order formulation:
$\partial_{x}\phi\simeq(\phi_{i+1,j}-\phi_{i-1,j})/2\delta x$ and
$\partial_{y}\phi\simeq(\phi_{i,j+1}-\phi_{i,j-1})/2\delta x$, where $i$ and
$j$ are the indexes of the coordinates $x$ and $y$, or the method based on the
directional derivatives with higher-order isotropy [54]. The isotropic finite-
difference approximations to the first-order derivative can be computed as:

$\nabla\phi=\frac{1}{c_{s}^{2}}\sum_{\alpha=0}^{Q}w_{\alpha}\left(|c_{\alpha}|^{2}\right)\phi\left(x+c_{\alpha}\right)c_{\alpha}$
(29)

where $w_{\alpha}\left(|c_{\alpha}|^{2}\right)$ are the weights associated to
each layer of neighboring nodes maximizing isotropy [55]. These weights are
summarized in table LABEL:t1. Given the importance of this directional
derivative in the growth dynamics of the solid phase, the choice of the
approximation will be briefly discussed here.

Table 1: Weights for 4th, 6th and 8th order isotropic tensors in two dimensions [56, 55] Tensor | $w(1)$ | $w(2)$ | $w(3)$ | $w(4)$ | $w(5)$ | $w(6)$ | $w(7)$ | $w(8)$
---|---|---|---|---|---|---|---|---
E4 | 1/3 | 1/12 | - | - | - | - | - | -
E6 | 4/15 | 1/10 | - | 1/120 | - | - | - | -
E8 | 4/21 | 4/45 | - | 1/60 | 2/315 | - | - | 1/5040
Figure 2: The grid points identifying the set of velocity fields for a 2d case
from [56, 55]. With reference to the weights reported in Table LABEL:t1,
different degrees of isotropy can be achieved: fourth order [up to
$\omega$(2)], sixth order [up to $\omega$(4)], eighth order [up to
$\omega$(8)].

To that end we consider a generic growing crystal with hexagonal symmetry
driven by the temperature field. For all simulations a domain of size
$800\times 800$, with $\delta_{x}=\delta_{y}=0.01$ and $\delta_{t}=1.5\times
10^{-5}$ is considered. Furthermore, $\tau_{0}=1.5625\times 10^{-4}$,
$W_{0}=0.0125$, $\kappa=1$, the undercooling is $\Delta=0.3$, $\lambda=10$,
$\varepsilon_{s}=0.05$. The circular seed is initialized at the center of the
domain $\bm{X_{c}}=(400,400)$ with a radius of $R_{s}=10$ lattice units.

Figure 3: Iso-contours of $\phi=0$ at time $t=1\times 10^{5}\delta_{t}$. Black
circles are reference results from [41].

The solid boundaries, $\phi=0$, as obtained using different approximations at
$t=10^{5}\delta_{t}$ are shown in Fig. 3. It could be seen that results
obtained with a second-order finite difference method (E2) do not match well
with the reference. As a function of their orientation relative to the main
axes, the dendrites are slightly different from each other. When the gradients
are calculated with higher order isotropic formulations these non-physical
anisotropic effects are considerably reduced. As such a minimum isotropy of
order 4, i.e. the E4 stencil, is necessary for a six-tip crystal simulation.

####  3.1.2 Validation for crystal growth

To validate the model and subsequent implementation, we use the generic system
with four-fold symmetry studied in [41, 53]. To provide both qualitative and
quantitative proof, we compare the shape of the dendrites and the evolution of
the tip velocity.
Initially, a circular seed of radius $R_{s}=10\delta_{x}$ is placed at the
centre of the square domain. The interface thickness $W_{0}$ and the
characteristic time $\tau_{0}$ are $W_{0}=\tau_{0}=1$. The grid-size is set to
$\delta_{x}/W_{0}=0.4$ [26] while the time-step size is $\delta_{t}=0.008$.
The capillary length $d_{0}$ and the kinetic coefficient $\beta$ are given by
[26]: $d_{0}=a_{1}W_{0}/\lambda$ and $\beta=a_{1}(\tau_{0}/\lambda
W_{0}-a_{2}W_{0}/\kappa)$ where $a_{1}=0.8839$ and $a_{2}=0.6267$. In this
benchmark, we choose the parameters $\kappa=4$,
$\lambda=\frac{\tau_{0}\kappa}{a_{2}W_{0}^{2}}=6.3826$ with $\beta=0$. The
anisotropic strength is $\epsilon_{s}=0.05$. The initial supersaturations are
$U_{0}=0.3$ and $0.55$. For $U_{0}=0.3$, we use a grid of size $1000^{2}$
nodes while for the latter the domain size is reduced to $500^{2}$ nodes. The
growth velocity of the crystal tips $\tilde{\bm{V_{p}}}$ is made dimensionless
as $\tilde{\bm{V_{p}}}=\bm{V_{p}}d_{0}/D$. The non-dimensional position
$\tilde{\bm{x}}$ is defined as ($\tilde{\bm{x}}=\bm{x}/\delta x$) and the
reduced time $T$ as ($T=t/\tau_{0}$).

Figure 4: (left) $\phi=0$ iso-contours for $U_{0}=0.3$ at $t=1.3\times
10^{5}\delta_{t}$. Red dots are from the present study while dashed black
lines are from [53]. (right) Dimensionless tip velocity $V_{p}$ as a function
of time (in units of $\tau_{0}$) for (lower curve/symbols) $U_{0}=0.3$ and
(upper curve/symbols) $U_{0}=0.55$. The red circles are from the present
study, while plain black lines are extracted from [53].

The obtained results are shown in Fig. 4. As shown there, the data obtained
from the present work closely follows those reported in [53].

####  3.1.3 Validation of flow/solid coupling

In order to validate the coupling between the flow field and other solvers, we
model the 2-D case presented in [43, 42], and solved there via an adaptive
finite-elements solver. The computational domain is a box of length $L=204.8$,
the grid-size $\delta_{x}=0.4$ and the time-step unit $\delta_{t}=0.008$. We
also set the interface thickness to $W_{0}=1$, the characteristic time to
$\tau_{0}=1$, the coupling strength to $\lambda=6.383$ and the anisotropy
strength to $\epsilon_{s}=0.05$. The kinematic viscosity of the flow field is
$\mu_{f}=92.4$. The flow enters from the left side of the box with a fixed
inlet velocity $u_{x}=1.0$. The boundary is set to zero-gradient. The top and
bottom sides of the box are set to periodic.
The evolution of all dendrite tips without/with velocity are plotted in Fig.
5.

Figure 5: Computed phase-field contours from the dendritic growth in 2-D
(left) without and (right) with flow at two different times. Red and blue
symbols are results at $t=72$ and $104$ while black dashed lines are
corresponding reference data.

It can be observed that the phase-field contours are in excellent agreement
with the reference. To better illustrate the interaction of the flow field
with the growing solid, the streamlines and species concentration field at two
different times are shown in Fig. 6. As expected, one can readily see the non-
isotropic distribution of the concentration field around the growing seed
caused by the flow field. The incoming velocity induces a higher concentration
gradient around the tip on the causing it to grow faster than its counter-part
in the opposite direction.

Figure 6: (top row) Velocity field streamlines and (bottom row) concentration
fields at two different times: (left) $T=72$ and (right) $104$.

###  3.2 Faceted and dendritic snowflake growth

In this subsection, we first study the plate growth regime of snowflakes as a
function of the supersaturation to showcase the ability of the proposed model
to capture the wide variety of habits exhibited by snowflakes. Then, we
briefly look at the effect of forced convection on the growth of snowflakes.

####  3.2.1 Ice crystal habit in thin-plate regime as a function of
temperature and supersaturation

As stated earlier, the focus of the present study is on the plate growth
regime. According to the Nakaya diagram, the widest variety of habits and
instabilities in this regime can be observed at $-16^{\circ}$C where the
saturated vapor density of ice is around
$\rho_{sat}^{I}=1.125\hbox{g}/\hbox{m}^{3}$ and the saturated vapor density of
water is around $\rho_{sat}^{W}=1.53\hbox{g}/\hbox{m}^{3}$ [53]. The excess
density over vapor-water equilibrium is $\Delta\rho=\rho-\rho_{sat}^{W}$ which
is exhibited in Fig. 1 by y-axis, where $\rho_{sat}^{W}$ is the satuaration
vapor density of water. The supersaturation is given by
$U=(c-c_{sat}^{I})/c_{sat}^{I}$, where $c_{sat}^{I}$ is the saturation number
density of vapor of ice at temperature $T$. Using $\rho=m_{H_{2}o}c$, where
$m_{H_{2}o}$ is the mass of a molecule of water, it could be written
$U=(\Delta\rho+\rho_{sat}^{W}-\rho_{sat}^{I})/\rho_{sat}^{I}$ and different
snowflake morphologies are obtained by changing the initial excess vapor
density $\Delta\rho$. According to [57], the melting temperature of the snow
crystals is $T_{m}=276.9\hbox{K}$, heat capacity $c_{p}=4.23\times
10^{6}\hbox{J}\hbox{m}^{-3}\hbox{K}^{-1}$, diffusion coefficient of the vapor
$D=1.17\times 10^{-7}\hbox{m}^{2}\hbox{s}^{-1}$ which is illusrated in Eq. 2,
surface tension $\gamma=2.845\times 10^{-2}\hbox{J}\hbox{m}^{-2}$, and latent
heat $L=1.12\times 10^{8}\hbox{J}\hbox{m}^{-3}$. The capillary length is
computed to be $d_{0}\simeq 2$ nm. The coupling parameter is set to
$\lambda=3$ [58].
The simulation parameters are chosen as: $W_{0}=1.25\delta_{x}$ and
$\tau_{0}=20\delta_{t}$. The simulations were conducted in a 400$\times$400
box. The anistropy strength was set to $\varepsilon_{xy}=0.05$. To produce
different snowflake morphologies, the following parameters were varied:
initial supersaturation $u_{0}$, depletion rate $L_{sat}$, initial density
$\rho_{0}$ and the excess density $\Delta\rho^{W}$. The corresponding values
are listed in table LABEL:t2.

Table 2: Parameters for different morphologies of snowflakes Nr. | shapes | $U_{0}[-]$ | $L_{sat}[-]$ | $\rho_{0}[\hbox{g}/\hbox{m}^{3}]$ | $\Delta\rho^{W}[\hbox{g}/m^{3}]$
---|---|---|---|---|---
1 | solid plate | 0.4 | 2.0 | 1.575 | 0.045
2 | stellar plate I | 0.5 | 1.8 | 1.6875 | 0.1575
3 | sectored plate | 0.5 | 1.6 | 1.6875 | 0.1575
4 | stars | 0.6 | 1.0 | 1.8 | 0.27
5 | fern dendrite | 0.8 | 1.6 | 2.025 | 0.495
6 | fernlike stellar | 0.7 | 1.6 | 1.9125 | 0.383
7 | stellar plate II | 0.5 | 1.2 | 1.6875 | 0.1575

The results obtained for the conditions listed in table LABEL:t2 are shown in
Fig. 7 along with corresponding numerical and experimental data from [27] and
[10].

Figure 7: The morphology numbers are from 1 to 7 (from left to right) in this
figure, as indicated in table LABEL:t2. Comparison between (top) real
snowflakes photographs taken from Libbrecht’s experiments [10], (middle) our
phase-field simulations, and (bottom) the simulations from [27].

The primary habit of the crystal (six-fold symmetry) is dictated by the
anisotropy function (and the microscopic crystallographic structure). At lower
supersaturation values where the adsorption rate is slow, the surface
diffusion process characteristic time is smaller and therefore dominates over
surface adsorption. More explicitly it means that the adsorbed water molecules
have enough time to propagate on the crystal surface and find the points with
the lowest potential (dictated by the molecules’ arrangement in the crystal
lattice). Furthermore, given the low growth rate and gradients, the surface is
not subject to branching instabilities. As the supersaturation goes up, the
larger adsorption rate at the sharper parts of the interface (regions with the
highest curvatures and consequently highest surface area) result in the
formation of six thick branches (usually referred to as primary branches). In
the lower supersaturation regimes these primary branches have a faceted
structure following the symmetry of the crystal. As the concentration goes
further up, the branches get thinner and rougher (the straight faces tend to
disappear); this eventually produces secondary instabilities and branches
going towards a somewhat fractal structure. All the obtained crystal habits,
are in excellent agreement with not only numerical simulations from [27] but
also experimental data from [10]. Further comparing the different crystal
habits to Nakaya’s diagram, it can be concluded that the proposed model
correctly predicts the behavior of the crystal in the platelet regime. The
next part will focus on the effects of ventilation on the evolution of the
crystal habit.

####  3.2.2 Ventilation effect on snowflakes

Growth of snowflakes under forced convection is a topic of interest, as for
example falling flakes are usually subject to ventilation. The dynamics of
snowflake growth under ventilation effects are not very well-documented.
Libbrecht and Arnold [59] worked on an aerodynamic model to show the growth
and appearance of new habits such as the triangular snow crystals both in
nature and in laboratory settings [60]. Furthermore, the anisotropy induced by
the flow field can, when pronounced, cause different regions on the surface of
the crystal to grow in different regimes. The present section will focus on
recovery of these effects and qualitative validation with experimental
observations.
For the first configuration, the domain has a size of $1600\times 1600$ grids,
the initial seed radius is $R=25\delta x$ and the initial vapor density
$\rho_{0}=1.364$g/m3 (supersaturation $U_{0}=0.2$). The grid size is set to
$\delta_{x}=4.8\times 10^{-6}$ m while the time-step is $\delta_{t}=2\times
10^{-7}$ s. The coupling strength is $\lambda=3$. The flow blows from the
bottom to the top along the $y$-axis. The inlet velocity is set to
$\bm{u}_{y}=0.12\hbox{m}/\hbox{s}$ and the kinetic viscosity to
$\nu_{f}=1.152\times 10^{-6}\hbox{m}^{2}/\hbox{s}$. The outlet is modeled
using a Neumann condition, while along the $x$-axis periodic boundary
conditions are used. The resulting snowflake morphology is shown in Fig. 8.

Figure 8: The morphology of the asymmetrical hexagonal snowflakes at $t=7500$
in units of $\tau_{0}$. Top left: experimental image, top right: our
simulation, and bottom: associated supersaturation field.

It can be clearly observed that the crystal growth on the side facing the
incoming flow is higher than its neighboring sides. Furthermore, its opposite
side is growing slower than its neighbour. Both of these non-symmetrical
growth rates push the habit towards, first a non-symmetrical hexagon and then
a triangular shape, in agreement with experimental observations.

Figure 9: (left) morphology of the snowflakes with ventilation effects and
(right) velocity field streamlines and supersaturation fields at (from top to
bottom) $t$ = 0, 4000 and 8000 in units of $\tau_{0}$.

To further put the effect of hydrodynamics into perspective we also consider
another test-case with the same configuration as the previous one, however
with a larger initial supersaturation. The initial vapor density is set to
$\rho_{f}=1.445\hbox{g}/\hbox{m}^{3}$ (supersaturation $U_{0}=0.3$). The inlet
velocity is set to $\bm{u}_{y}=0.24\hbox{m}/\hbox{s}$. The evolution of the
crystal, stream-lines and supersaturation fields at different time-steps are
illustrated in Fig. 9.
As observed in these figures, the natural anisotropy in supersaturation around
the crystal is further accentuated by the formation of two re-circulation
zones. Due to the presence of these flow structures, the growth rate on the
top half of the crystal is slowed down and brought into the regular hexagonal
habit regime. The lower half, however, is subject to larger concentration
gradients (further amplified by the incoming convective flux) and therefore
exhibits primary branching instabilities. As the system progresses further in
time, the bottom-facing branches distinguish themselves from the side-
branches, as we observe secondary instability effects on them, along with a
much faster growth rate on the main branches. These observations are in clear
agreement with expectations from fluid dynamics. Going into the details of the
crystal structure and supersaturation fields, we can also see that secondary
branching instabilities are not present at the bases of the down-facing
primary branches. That is explained by the fact that they are in a flow
stagnation zone (due to the built-up pressure in this closed area), where the
supersaturation is almost fully depleted. All of these effects are in
qualitative agreement with experimental observations reported in the
literature [34].

##  4 Conclusions and perspectives

In this study, we presented a model for the simulation of snowflakes.
Throughout the many generic test-cases it has been shown that this model
closely matches results from other solvers for the system of macroscopic s
needed for this system. It was further shown that the proposed formulation is
able to capture the different crystal habits in the plate growth regime. The
final crystal habits were in good agreement with experimental data and the
Nakaya diagram predictions [9]. To go one step further, the model was also
used to look at the possible effects of forced convection on the growth
dynamics and resulting asymmetrical shapes.
Given the promising results obtained using this model, future work will focus
on the extension of this study to the long-prismatic growth regime to cover
the entire spectrum of habits exhibited by snowflakes. Effects from local
variations in temperature will also be added to the model to have a better
image of the mechanisms behind the growth of snowflakes [61].

##  5 Acknowledgement

Q.T. would like to acknowledge the financial support by the EU-program ERDF
(European Regional Development Fund) within the Research Center for Dynamic
Systems (CDS). S.A.H. would like to acknowledge the financial support of the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in TRR 287
(Project-ID 422037413), and thank K. G. Libbrecht for allowing the authors to
use images from his ice crystal growth experiments.

## References

  * [1] F. Frank, Snow crystals, Contemporary Physics 23 (1) (1982) 3–22.
  * [2] B. J. Mason, Snow crystals, natural and man made, Contemporary Physics 33 (4) (1992) 227–243.
  * [3] P. K. Wang, Shape and microdynamics of ice particles and their effects in cirrus clouds, Advances in Geophysics 45 (2002) 1–258.
  * [4] I. Shafranovskii, Kepler’s Crystallographic ideas and his tract “the six-cornered snowflake”, Vistas in Astronomy 18 (1975) 861–876.
  * [5] F. E. Chickering, Cloud Crystals: A Snow-flake Album, Appleton, 1864.
  * [6] W. Bentley, W. Humphreys, Snow crystals, McGraw-Hili Book Company, Journal of Interdisciplinary History.
  * [7] U. Nakaya, Snow crystals: Natural and Artificial, Harvard University Press, 2013.
  * [8] T. Gonda, S. Nakahara, Dendritic ice crystals with faceted tip growing from the vapor phase, Journal of Crystal Growth 173 (1-2) (1997) 189–193.
  * [9] U. Nakaya, The formation of ice crystals, Compendium of Meteorology (1951) 207–220.
  * [10] K. G. Libbrecht, Snowcrystals.com (1999).
URL http://www.snowcrystal.com

  * [11] U. Nakaya, Snow crystals, natural and artificial, Harvard University Press (1954) 510.
  * [12] K. G. Libbrecht, The physics of snow crystals, Reports on Progress in Physics 68 (4) (2005) 855.
  * [13] K. G. Libbrecht, Morphogenesis on ice: The physics of snow crystals, Engineering and Science 64 (1) (2001) 10–19.
  * [14] K. G. Libbrecht, Physical dynamics of ice crystal growth, Annual Review of Materials Research 47 (2017) 271–295.
  * [15] G. Ross, The snowflake: Winter’s Secret Beauty, American Scientist 92 (2) (2004) 181–182.
  * [16] B. Mason, The growth of ice crystals from the vapour and the melt, Advances in Physics 7 (26) (1958) 235–253.
  * [17] J. Nelson, C. Knight, Snow crystal habit changes explained by layer nucleation, Journal of the Atmospheric Sciences 55 (8) (1998) 1452–1465.
  * [18] J. Nelson, Growth mechanisms to explain the primary and secondary habits of snow crystals, Philosophical Magazine A 81 (10) (2001) 2337–2373.
  * [19] K. G. Libbrecht, M. E. Rickerby, Measurements of surface attachment kinetics for faceted ice crystal growth, Journal of Crystal Growth 377 (2013) 1–8.
  * [20] C. Magono, C. Lee, Meteorological classification of natural snow crystals, Journal of the Faculty of Science 2 (4) (1966) 321–335.
  * [21] J. Gravner, D. Griffeath, Modeling snow-crystal growth: A three-dimensional mesoscopic approach, Physical Review E 79 (1) (2009) 011601.
  * [22] J. G. Kelly, E. C. Boyer, Physical improvements to a mesoscopic cellular automaton model for three-dimensional snow crystal growth, Crystal Growth & Design 14 (3) (2014) 1392–1405.
  * [23] D. Chatterjee, S. Chakraborty, A hybrid lattice boltzmann model for solid–liquid phase transition in presence of fluid flow, Physics Letters A 351 (4-5) (2006) 359–367.
  * [24] J. W. Barrett, H. Garcke, R. Nürnberg, Numerical computations of faceted pattern formation in snow crystal growth, Physical Review E 86 (1) (2012) 011604\.
  * [25] I. Singer-Loginova, H. M. Singer, The phase field technique for modeling multiphase materials, Reports on Progress in Physics 71 (10) (2008) 106501\.
  * [26] A. Karma, W. J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Physical Review E 57 (4) (1998) 4323.
  * [27] G. Demange, H. Zapolsky, R. Patte, M. Brunel, A phase field model for snow crystal growth in three dimensions, NPJ Computational Materials 3 (1) (2017) 1–7.
  * [28] G. Demange, H. Zapolsky, R. Patte, M. Brunel, Growth kinetics and morphology of snowflakes in supersaturated atmosphere using a three-dimensional phase-field model, Physical Review E 96 (2) (2017) 022803.
  * [29] G. D. Doolen, Lattice boltzmann methods for fluid flows, Annual Review of Fluid Mechanics 30 (1) (1998) 329–364.
  * [30] H. W. Zheng, C. Shu, Y. T. Chew, A lattice boltzmann model for multiphase flows with large density ratio, Journal of Computational Physics 218 (1) (2006) 353–371.
  * [31] S. A. Hosseini, N. Darabiha, D. Thévenin, Mass-conserving advection–diffusion lattice boltzmann model for multi-species reacting flows, Physica A: Statistical Mechanics and its Applications 499 (2018) 40–57.
  * [32] S. A. Hosseini, H. Safari, N. Darabiha, D. Thévenin, M. Krafczyk, Hybrid lattice boltzmann-finite difference model for low mach number combustion simulation, Combustion and Flame 209 (2019) 394–404.
  * [33] S. A. Hosseini, A. Abdelsamie, N. Darabiha, D. Thévenin, Low-mach hybrid lattice boltzmann-finite difference solver for combustion in complex flows, Physics of Fluids 32 (7) (2020) 077105.
  * [34] D. Medvedev, K. Kassner, Lattice boltzmann scheme for crystal growth in external flows, Physical Review E 72 (5) (2005) 056703.
  * [35] I. Rasin, W. Miller, S. Succi, Phase-field lattice kinetic scheme for the numerical simulation of dendritic growth, Physical Review E 72 (6) (2005) 066705\.
  * [36] W. Miller, I. Rasin, S. Succi, Lattice boltzmann phase-field modelling of binary-alloy solidification, Physica A: Statistical Mechanics and its Applications 362 (1) (2006) 78–83.
  * [37] D. Medvedev, T. Fischaleck, K. Kassner, Influence of external flows on crystal growth: Numerical investigation, Physical Review E 74 (3) (2006) 031606.
  * [38] C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice boltzmann model for melting with natural convection, International Journal of Heat and Fluid Flow 29 (5) (2008) 1469–1480.
  * [39] D. Sun, M. Zhu, S. Pan, D. Raabe, Lattice boltzmann modeling of dendritic growth in a forced melt convection, ACTA MATERIALIA 57 (6) (2009) 1755–1767.
  * [40] G. Lin, J. Bao, Z. Xu, A three-dimensional phase field model coupled with a lattice kinetics solver for modeling crystal growth in furnaces with accelerated crucible rotation and traveling magnetic field, Computers & Fluids 103 (2014) 204–214.
  * [41] A. Younsi, A. Cartalade, On anisotropy function in crystal growth simulations using lattice boltzmann equation, Journal of Computational Physics 325 (2016) 1–21.
  * [42] J. H. Jeong, N. Goldenfeld, J. A. Dantzig, Phase field model for three-dimensional dendritic growth with fluid flow, Physical Review E 64 (4) (2001) 041602.
  * [43] C. Beckermann, H. J. Diepers, I. Steinbach, A. Karma, X. Tong, Modeling melt convection in phase-field simulations of solidification, Journal of Computational Physics 154 (2) (1999) 468–496.
  * [44] B. Echebarria, R. Folch, A. Karma, M. Plapp, Quantitative phase-field model of alloy solidification, Physical Review E 70 (6) (2004) 061604.
  * [45] R. Folch, M. Plapp, Publisher’s note: Quantitative phase-field modeling of two-phase growth [phys. rev. e 72, 011602 (2005)], Physical Review E 72 (2) (2005) 029903.
  * [46] S. A. Hosseini, N. Darabiha, D. Thévenin, Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature, Philosophical Transactions of the Royal Society A 378 (2175) (2020) 20190399.
  * [47] X. Shan, X. F. Yuan, H. Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation, Journal of Fluid Mechanics 550 (2006) 413–441.
  * [48] S. A. Hosseini, C. Coreixas, N. Darabiha, D. Thévenin, Extensive analysis of the lattice boltzmann method on shifted stencils, Physical Review E 100 (6) (2019) 063301.
  * [49] S. A. Hosseini, Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows, Ph.D. thesis, Université Paris-Saclay & Otto-von-Guericke university (2020).
  * [50] S. Ponce D., S. Chen, G. D. Doolen, Lattice boltzmann computations for reaction-diffusion equations, The Journal of Chemical Physics 98 (2) (1993) 1514–1523.
  * [51] C. Huber, B. Chopard, M. Manga, A lattice boltzmann model for coupled diffusion, Journal of Computational Physics 229 (20) (2010) 7956–7976.
  * [52] S. D. Walsh, M. O. Saar, Macroscale lattice-boltzmann methods for low peclet number solute and heat transport in heterogeneous porous media, Water Resources Research 46 (7).
  * [53] A. Cartalade, A. Younsi, M. Plapp, Lattice boltzmann simulations of 3d crystal growth: Numerical schemes for a phase-field model with anti-trapping current, Computers & Mathematics with Applications 71 (9) (2016) 1784–1798.
  * [54] L. Chen, Q. Kang, Y. Mu, Y. He, W. Tao, A critical review of the pseudopotential multiphase lattice boltzmann model: Methods and applications, International Journal of Heat and Mass Transfer 76 (2014) 210–236.
  * [55] X. Shan, Analysis and reduction of the spurious current in a class of multiphase lattice boltzmann models, Physical Review E 73 (4) (2006) 047701.
  * [56] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, Generalized lattice boltzmann method with multirange pseudopotential, Physical Review E 75 (2) (2007) 026702.
  * [57] I. Yoshizaki, T. Ishikawa, S. Adachi, E. Yokoyama, Y. Furukawa, Precise measurements of dendrite growth of ice crystals in microgravity, Microgravity Science and Technology 24 (4) (2012) 245–253.
  * [58] J. C. Ramirez, C. Beckermann, A. s. Karma, H. J. Diepers, Phase-field modeling of binary alloy solidification with coupled heat and solute diffusion, Physical Review E 69 (5) (2004) 051607.
  * [59] K. G. Libbrecht, Physically derived rulesfor simulating faceted crystal growth using cellular automata, ArXiv Condensed Matter: Materials Science.
  * [60] K. G. Libbrecht, H. M. Arnold, Aerodynamic stability and the growth of triangular snow crystals, ArXiv Condensed Matter: Materials Science.
  * [61] S. A. Hosseini, N. Darabiha, D. Thévenin, Lattice boltzmann advection-diffusion model for conjugate heat transfer in heterogeneous media, International Journal of Heat and Mass Transfer 132 (2019) 906–919.