2022

[1]Denys Y. Kononenko [1]Oleg Janson

1]Institute for Theoretical Solid State Physics, Dresden, 01069, Germany
2]Institute for Theoretical Physics, TU Dresden, Dresden, 01069, Germany

# Data-driven estimation of transfer integrals in undoped cuprates

<EMAIL_ADDRESS>Ulrich K. Rößler  Jeroen van den Brink
<EMAIL_ADDRESS>[  [

###### Abstract

Undoped cuprates are an abundant class of magnetic insulators, in which the
synergy of rich chemistry and sizable quantum fluctuations leads to a variety
of magnetic behaviors. Understanding the magnetism of these materials is
impossible without the knowledge of the underlying spin model. The typically
dominant antiferromagnetic superexchanges can be accurately estimated from the
respective electronic transfer integrals. Density functional theory
calculations mapped onto an effective one-orbital model in the Wannier basis
are an accurate, albeit computationally cumbersome method to estimate such
transfer integrals in cuprates. We demonstrate that instead an Artificial
Neural Network (ANN), trained on the results of high-throughput calculations,
can predict the transfer integrals using the crystal structure as the only
input. Descriptors of the ANN model encode the spatial configuration and the
chemical composition of the local crystalline environment. A virtual toolbox
employing our model can be readily employed to determine leading superexchange
paths as well as for rapidly assessing the relevant spin model in yet unknown
cuprates.

###### keywords:

quantum magnetism, machine learning, high-throughput calculations, transfer
integrals

##  1 Introduction

The data-driven approach accompanied by modern machine-learning (ML)
techniques becomes an increasingly important tool of scientific investigations
across many domains of physics. From quantum to fluid mechanics tsubaki.prl.20
; brunton.annurev-fluid.20  learning from data facilitates descriptions of
complex phenomena for which analytical approaches are prohibitively
challenging. The adoption of ML in solid-state physics and material science is
particularly appealing: the sheer amount of collected experimental and
computed records propels the community to design data-driven frameworks for
prediction of materials properties scmidt.npj.compmat.19 ; xie.prl.20 .

The application of ML for problems of solid-state physics is not
straightforward. One of the key challenges is to represent periodic
(crystalline) or finite (molecule, local crystal environment, etc.) atomic
systems as descriptors – data structures amenable to ML methods. Such
descriptors must be invariant with respect to the choice of the unit cell
(crystals) or to global rotations in a finite system. Several classes of
descriptors have been developed for material properties prediction: Coulomb
matrix rupp.prl.12 , partial radial distribution function (PRDF)
schtt.prb.2014 , smooth overlap of atomic positions (SOAP) bartok.prb.2013 ,
diffraction fingerprint (DF) ziletti.natcomm.2018  and three-dimensional (3D)
Zernike descriptor (3DZRD) venkatraman.joc.09 . The latter were designed for
the characterization of 3D shapes novotni2003 ; novotni.shape.04  and
successfully employed for comparison of molecules sael.prot.2008 ;
mak.jmgm.2008 ; sael2008 . These descriptors are invariant with respect to the
number of chemical species in the dataset, they store detailed information
about the spatial arrangement, and do not require additional simulation
software. These features as well as the compact size of the resulting data
structures make 3DZRD ideally suited for a description of diverse, dissimilar
crystalline environments.

Another key challenge is the construction of descriptors that represent
material properties. For instance, it is widely accepted that electronic,
magnetic, and topological properties of bulk materials are rooted – in a
highly nontrivial way — in their electronic structure. However, using the
complete band structure of a material as a universal descriptor is not
possible for a number of reasons: different number of bands, non-universal
discretization of the Brillouin zone, huge dimensionality etc. A more
practical approach is to restrict the description to the states relevant for
the physical quantity of interest. Naturally, this is possible only for a
certain class of materials and only for specific physical property.

Following this idea, we apply a data-driven approach to assess spin models in
undoped cuprates – stoichiometric inorganic materials containing divalent
copper and oxygen atoms. In contrast to their doped counterparts – the high-
temperature cuprate superconductors Plakida2010htc  – undoped cuprates are
magnetic insulators with the $3d^{9}$ electronic configuration of Cu2+. The
sizable Jahn-Teller distortion lifts the orbital degeneracy, giving rise to
half filling and localized $S$=$\frac{1}{2}$ spins. Owing to the plethora of
structure types and the quantum limit assured by $S$=$\frac{1}{2}$, undoped
cuprates exhibit a variety of magnetic behaviors sahadasgupta2021 , from
simple quantum dimers and spin chains – to exotic collective behaviors such as
the spin-liquid regime in herbertsmithtite $\gamma$-Cu3Zn(OH)6Cl2 helton07 ,
Bose-Einstein condensation of magnons in Han purple BaCuSi2O6 jaime04 , or
bound magnon states in volborthite Cu3V2O7(OH)${}_{2}\cdot 2$H2O kohama19 .

Understanding the magnetic properties of cuprates requires the knowledge of
the underlying spin model. While exchange anisotropies are generally present
and can alter the magnetic properties, the backbone of spin models are
isotropic interactions, and the relevant Heisenberg Hamiltonian is the
following sum:

$\displaystyle\mathcal{H}=\frac{1}{2}\sum_{ij}J_{ij}\left(\mathbf{S}_{i}\cdot\mathbf{S}_{j}\right),$
(1)

where $\mathbf{S}_{i}$ and $\mathbf{S}_{j}$ are spin operators on sites $i$
and $j$. The set of relevant magnetic exchange integrals $\\{J_{ij}\\}$
determines the spin model. It is important to note that each individual
$J_{ij}$ term is a sum of antiferromagnetic ($J^{\text{AF}}_{ij}<0$) and
ferromagnetic ($J^{\text{FM}}_{ij}<0$) contributions that are driven by
competing processes Goodenough_Magnetism . Commonly,
$J^{\text{AF}}_{ij}\gg\lvert{}J^{\text{FM}}_{ij}\rvert$, with the exception of
short-range exchanges for which the ferromagnetic contribution can become
dominant.

The antiferromagnetic contribution is a textbook example of the superexchange
mechanism and can be derived via second-order perturbation theory of the
Hubbard model in the strong-coupling limit at half-filling as
$J^{\textsc{AF}}_{ij}=4t_{ij}^{2}/U_{\text{eff}}$ hubbard63 ; anderson59 ;
auerbach2000 . Here, $U_{\text{eff}}$ is the Coulomb repulsion within an
effective molecularlike orbital, which in most cuprates is dominated by
$3d_{x^{2}-y^{2}}$ orbital of Cu and $\sigma$-bonded $2p$ orbitals of O. There
is empirical evidence that $U_{\text{eff}}$ from the range 4–5 eV gives a
proper description of the magnetism of cuprates belik2004 ; johannes2006 ;
janson2009 . Hence, the knowledge of transfer integrals $t_{ij}$ paves the way
to a quantitative assessment of the spin model in the majority of cuprate
materials. Yet, extracting $t_{ij}$ directly from the structural information
is essentially impossible; instead, it requires first-principles calculations
followed by an additional modeling.

To overcome this challenge, we propose a data-driven approach for prediction
of transfer integrals in cuprates, which requires the crystal structure as the
only input. Our approach is based on the local crystal environment description
utilizing 3DZRD. The local crystal environment descriptor is used as input for
the ML model which is trained on the results of high-throughput density-
functional-theory (DFT) calculations for hundreds of cuprate materials. DFT
calculation for each material is followed by automatized Wannierization and a
manual quality control. The trained model is wrapped into a freely accessible
web application 111https://smc-t.ifw-dresden.de/ that can be used for a quick
estimation of relevant transfer paths in new cuprate materials.

The paper is organized as follows: Section 2 describes the high-throughput DFT
calculations and the dataset of transfer integrals. Section 3 details the
descriptor for the Cu..Cu bonds. In Section 4 we compare three different ML
approaches for predicting transfer integrals and estimate the accuracy by a
cross-validation procedure (CV). In the last section, we discuss the
performance of our ML model for different classes of cuprates. In particular,
for the parent (undoped) compounds of high-temperature superconducting
cuprates we show that the ANN model quantitatively captures the ratio between
nearest and next-nearest-neighbor transfer integrals.

##  2 Dataset generation

We start with the description of the high-throughput DFT calculations employed
for the generation of the dataset of transfer integrals. The list of materials
contains 672 unique structures of undoped cuprates. The structures were
filtered out from the 10 710 cuprate structures stored in the Inorganic
Crystal Structure Database (ICSD) bergerhoff.crystallographic.1987 . For this
screening, the following criteria were consecutively applied: (i) the presence
of Cu2+ ions, (ii) electroneutrality (zero total charge), (iii) absence of
sites with fractional occupancies, (iv) the minimal inter-atomic distance of
0.5 Å, and (v) the absence of other magnetic atoms beyond Cu suppl . The
latter criterion is necessary to filter out compounds with multiple magnetic
atoms, where the presence of additional bands in the relevant energy range may
render the effective one-orbital model inapplicable and its results
misleading. For the analysis of the crystal structures we used the pymatgen
library Ong2013pymatgen  for Python.

For each structure, we performed DFT calculations to construct Wannier
Hamiltonians and subsequently determined the transfer integrals. All DFT
calculations were performed using the generalized gradient approximation (GGA)
perdew.prl.1996  with the full potential code FPLO of version 18.00-52
koepernik.prb.1999 . The computational workflow comprised several steps.
First, scalar-relativistic nonmagnetic DFT calculations were carried out and
the Hellmann-Feynman forces were calculated. Second, for hydrogen-containing
compounds whose calculated forces exceeded the threshold of 0.1 eV/Å, we
optimized the internal coordinates of H atoms within GGA. The rationale behind
this step are largely inaccurate H positions as determined by x-ray
diffraction (which is by far most common method of structure determination).
Since a considerable number of cuprates contain hydrogen, typically as
hydroxyl groups or water molecules, inclusion of such partly optimized
structures allowed us to considerably extend the data set. All other cuprates
whose forces exceeded the threshold 0.1 eV/Å were discarded. Third, we
calculated the orbital-resolved density of states (DOS) and band structure.
From the orbital-resolved DOS, the energy interval which contains the copper
$3d_{x^{2}-y^{2}}$ bands was determined. The energy interval is selected such
that the contribution of the magnetic $3d_{x^{2}-y^{2}}$ orbital in the total
density of Cu $3d$ states exceeds 5 %. The determined energy window
$[\mathscr{E}_{\textsc{min}},\mathscr{E}_{\textsc{max}}]$ was adopted for
Wannierization in the next step. Fourth, the Wannier transformation procedure
was performed to obtain the effective one-orbital Hamiltonian $H$ in the
Wannier basis. We used copper $3d_{x^{2}-y^{2}}$ orbitals as projectors and
the interval $[\mathscr{E}_{\textsc{min}},\mathscr{E}_{\textsc{max}}]$ as the
energy window to construct the Wannier functions (WF). The latter is necessary
to discriminate the target antibonding orbital (crossing the Fermi energy)
from its bonding sibling at the bottom of the valence band. The transfer
integral between two WF $w_{i}$ and $w_{j}$ placed at copper sites $i$, $j$ is
determined as (real) Hamiltonian matrix element
$t_{ij}=\langle{}w_{i}\lvert{}H\rvert{}w_{j}\rangle$. The details on the
construction of Wannier functions and the construction of the respective
tight-binding models are provided in the papers koepernik2021arxiv ;
eschrig2009 .

After the calculation pipeline was completed, we obtained a list of transfer
integrals $\\{t_{ij}\\}$ that connected $i$-th and $j$-th copper sites
situated at the distance $r_{ij}$ from each other for all valid structures
suppl . For construction of the dataset we select transfer integrals larger
than 5 meV with Cu..Cu spacing less than 8 Å. The distribution of calculated
transfer integrals $t_{ij}$ among Cu..Cu distances is shown in Fig. 2. The
crystal chemistry of cuprates sets a natural lower limit for the bond lengths;
accordingly, there is no transfer integral with the distance less than 2.4 Å
in the dataset. Remarkably, for the vast majority of transfer integrals, the
absolute values are below 0.2 eV. This natural imbalance of the dataset will
inevitably affect the performance of predictions.

Figure 1: (Color online) Transfer integrals obtained from the DFT calculations
as a function of the Cu..Cu distance. Region I harbors transfer integrals
between edge-sharing (a) and corner-sharing (b) CuO4 plaquettes, while region
II is dominated by transfer integrals between CuO4 plaquettes that do not
share oxygen atoms (c). The inset shows the distribution of the computed
transfer integrals.

Figure 2: (Color online) Schematics of the workflow: selection of the local
crystal environment from the cuprate crystal structure, generation of
rotationally invariant descriptor $\vec{D}$ via decomposition of the local
crystal environment function in the truncated basis of 3DZF and prediction of
the transfer integral $t_{ij}$ with ML algorithm trained on the dataset from
DFT calculations. The illustrating example Ba2CuHgO4 (ICSD Identifier 75724)
hosts pairs of corner-sharing CuO4 square-like plaquettes.

##  3 Crystal Environment Descriptor

To describe the crystal environment, we first determine the midpoint $\vec{p}$
between a given pair of copper atoms and build a sphere with the empirically
determined threshold radius $R^{\textsc{th}}=\max(4,r_{ij}/2+0.2)$ Å centered
at $\vec{p}$. Next, all atoms in the sphere are enlisted in the crystal
environment alongside with nearest neighbors of $i$ and $j$. We consider
nearest neighbors as atoms distanced from $i$ or $j$ not farther than 2.5 Å.
After the local crystal environment is assembled, we shift the coordinate
system origin to the centroid (the point between copper pair) $\vec{p}$ and
normalize atoms coordinates by $r_{0}=6$ Å to fit the crystal environment in
the unit ball. To construct a robust representation of the local crystal
environment we introduce the piecewise function of site positions
$\mathcal{I}(\vec{r})$. The function $\mathcal{I}$ equal to the $q$-th atom
oxidation number $O_{q}$ in the ball with center at the position of $q$-th
atom $\vec{r}_{q}$ and radius $R_{q}$ equals to the ionic radius of the atom

$\displaystyle\mathcal{I}(\vec{r})=\begin{cases}O_{q}&\|\vec{r}_{q}-\vec{r}\|\leq
R_{q},\\\ 0&\text{otherwise}.\end{cases}$ (2)

The normalization factor $r_{0}$ is a sum of the maximal considered Cu-Cu
distance $\max\|\vec{r}_{ij}\|=4$ Å and a maximal considered ionic radius $2$
Å. The function $\mathcal{I}(x,y,z)$ describes the spatial configuration and
chemical composition of the crystal environment placed in the unit ball with
$x^{2}+y^{2}+z^{2}\leq 1$. An example of the local crystal environment defined
by (2) is shown in Fig. 2.

We describe the selected crystal environment $\mathcal{I}$ in the form of a
finite-dimensional vector. Such representation provides a robust way for
numerical operations with crystal environments, e.g. similarity and sorting.
To obtain the finite vector representation of the crystal environment we
decompose the $\mathcal{I}(x,y,z)$ in the truncated basis of three-dimensional
(3D) Zernike functions (3DZF) $Z^{m}_{nl}$ which are defined as follows
canterakis96 ; canterakis99 ; novotni.shape.04

$\displaystyle\begin{split}Z_{nl}^{m}(r,\theta,\phi)&=R_{nl}(r)Y_{lm}(\theta,\phi),\\\
R_{nl}(r)&=\sum_{\nu=0}^{(n-l)/2}Q_{l\nu}r^{2\nu+l},\\\
Q_{l\nu}&=\dfrac{(-1)^{k+\nu}}{4^{k}}\sqrt{\dfrac{2l+4k+3}{3}}\dfrac{\binom{2k}{k}\binom{k}{\nu}\binom{2(k+l+\nu)+1}{2k}}{\binom{k+l+\nu}{k}},\end{split}$
(3)

where indices $n$ and $l$ are positive integers which satisfy condition $n\geq
l$; $m$ changes from $-l$ to $l$ with constraint $(n-l)$ is even number;
$k=(n-l)/2$ and $Y_{lm}(\theta,\phi)$ are spherical harmonics, and
$(r,\theta,\phi)$ are spherical coordinates nist.math.handbook . For
convenience, we use Cartesian coordinates $(x,y,z)$ representation of 3DZF
implying change of coordinates:
$\mathcal{Z}_{nl}^{m}(x,y,z)=Z_{nl}^{m}\left(\sqrt{x^{2}+y^{2}+z^{2}},\arctan{\sqrt{x^{2}+y^{2}}/z},\arctan{y/x}\right)$.
3DZF form the complete basis of orthogonal functions in the unit ball, so that
the function $\mathcal{I}(x,y,z)$ defined in the unit ball
$x^{2}+y^{2}+z^{2}\leq 1$ can be expanded in the introduced basis
morais.mathematics-of-computation.2014 .

The decomposition coefficients read

$\displaystyle
c_{nl}^{m}=\dfrac{1}{V}\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}\mathcal{I}(x,y,z)\mathcal{Z}_{nl}^{m}\left(x,y,z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z,$
(4)

where the normalization factor is the volume of the unit ball $V=4\pi/3$ suppl
.

Note, $c_{nl}^{m}$ is not invariant with respect to rotations of the crystal
environment $\mathcal{I}$. Rotationally invariant characteristics can be
obtained by assembling the vector $\vec{C}_{nl}$ whose components are all
$(2l+1)$ coefficients with different $m$ for given pair of $n$ and $l$. The
norm of obtained vector $\|\vec{C}_{nl}\|=C_{nl}$ determined as

$\displaystyle C_{nl}=\left\|c_{nl}^{0},...,c_{nl}^{l-1},c_{nl}^{l}\right\|,$
(5)

is invariant with respect to the rotation of the crystal environment, thus the
pre-alignment is not required.

We introduce the finite dimensional vector-descriptor of the crystal
environment $\mathcal{I}$ as

$\displaystyle\vec{D}=\left(C_{00},C_{11},C_{20},...,C_{n_{\textsc{max}}l_{\textsc{max}}},r_{ij}\right),$
(6)

where copper-copper distance $r_{ij}$ is incorporated into the descriptor as
well. The size of the descriptor $\vec{D}$ is determined by the cut-off order
of the Zernike 3D moments $n_{\textsc{max}}$ and corresponding
$l_{\textsc{max}}$ in the truncated basis. The size of the 3DZF basis grows
with $n_{\textsc{max}}$ as the sum of the series
$\sum_{n=0}^{n_{\textsc{max}}}(n^{2}+3n+2)/2$. In the present work, we chose
the cut-off order $n_{\textsc{max}}=25$. The vector $\vec{D}$ encodes the
information about spatial configuration and chemical composition of the
crystal environment, allowing the introduction of the mapping of $\vec{D}$ on
the transfer integral.

Figure 3: (Color online) Performance of the ensemble ANN model on the testing
and training datasets for random split. The ensemble ANN shows $R^{2}=0.7$,
RMSE = 28 meV and MAE = 18 meV on the testing dataset and $R^{2}=0.9$, RMSE =
19 meV and MAE = 11 meV on the training dataset. The inset figure shows the
distribution of the ensemble ANN model error for the test dataset with $\mu$
and $\sigma$ are mean value and standard deviation of the errors. The solid
line corresponds to the normal distribution with parameters $\mu$ and
$\sigma$.

##  4 Transfer Integral Prediction

Our high-throughput DFT calculations yielded $N=1800$ local crystal
environments $\\{\vec{D}\\}$ with corresponding transfer integrals
$\\{t_{ij}\\}$. We build the ML model to predict the continuous-valued
attribute $t_{ij}$ associated with the local crystal environment descriptor
$\vec{D}$. To solve this regression problem, we tested the following models:
(i) linear (LIR), (ii) random forest (RFR) Breiman_2001  regression models,
and (iii) ANN. To achieve robust generalization and stability of the models,
we employ the bagging (bootstrap aggregating) ensemble technique
breiman1996bagging . The main idea behind bagging is to train multiple
instances of the same model on different subsets of the training data and then
combine their predictions to make the final estimation. For each model in the
ensemble, a random sample is drawn with a replacement from the original
training dataset. Thus, some data points may appear multiple times in the
sample, while others may be left out. When making predictions, the individual
predictions from each model are combined using a voting method (for
classification tasks) or averaging (for regression tasks). In the work, we
employ ensemble models with 100 estimators.

As a metric for the ensemble regression model performance with predictions
$\tau$ we use: (i) the coefficient of determination

$\displaystyle R^{2}=1-\dfrac{S_{\textsc{reg}}}{S_{\textsc{tot}}},$ (7)

where $S_{\textsc{reg}}=\sum_{p=1}^{M}(t^{p}_{ij}-\tau^{p})^{2}$ is a sum of
squared residuals of the regression model and
$S_{\textsc{tot}}=\sum_{p=1}^{M}(t^{p}_{ij}-\overline{t_{ij}})^{2}$ is a total
sum of squares with $\overline{t_{ij}}$ being the mean value of transfer
integral in the test dataset with $M$ samples.

(ii) the root mean squared error (RMSE)

$\displaystyle\text{RMSE}=\sqrt{\dfrac{1}{M}\sum_{p=1}^{M}(t^{p}_{ij}-\tau^{p})^{2}},$
(8)

and (iii) mean absolute error

$\text{MAE}=\dfrac{1}{M}\sum_{p=1}^{M}|t^{p}_{ij}-\tau^{p}|.$ (9)

For RMSE, squared errors of the model are included in the average, making this
measure more sensitive to outliers. As more variance in predictions, a larger
RMSE. The MAE provides a mean of linear scores with all errors weighted
equally.

Model | $\overline{R^{2}}$ |  $\overline{\text{RMSE}}$, meV |  $\overline{\text{MAE}}$, meV
---|---|---|---
LIR | 0.27 $\pm$ 0.18 | 44 $\pm$ 5 | 28 $\pm$ 1
RFR | 0.59 $\pm$ 0.01 | 34 $\pm$ 3 | 21 $\pm$ 1
ANN | 0.69 $\pm$ 0.05 | 29 $\pm$ 3 | 18 $\pm$ 1
Table 1: Results of shuffle-split CV of the selected ensemble models. The average value of $R^{2}$, RMSE, and MAE on six splits is given alongside the standard deviation. Model | $\overline{R^{2}}$ |  $\overline{\text{RMSE}}$, meV |  $\overline{\text{MAE}}$, meV
---|---|---|---
LIR | 0.27 $\pm$ 0.20 | 43 $\pm$ 6 | 28 $\pm$ 1
RFR | 0.57 $\pm$ 0.06 | 34 $\pm$ 5 | 21 $\pm$ 1
ANN | 0.69 $\pm$ 0.10 | 28 $\pm$ 3 | 18 $\pm$ 1
Table 2: Results of $k$-fold CV of the selected ensemble models with six
folds. The average value of $R^{2}$, RMSE, and MAE is given alongside the
standard deviation.

For model selection, we employ two CV strategies: the shuffle-split and
$k$-fold. In the shuffle-split, the dataset is randomly shuffled and then
split into training and test subsets containing a specific percentage of the
original data. The procedure is repeated the specified number of iterations.
In each iteration the model is trained and evaluated accordingly. The shuffle-
split procedure was implemented using the scikit-learn library scikit-learn
with six splits, a test size of 20 %, and a training size of 80 % of the
entire dataset. The results of the CV are presented in the Table 1. The model
selection procedure shows that ensemble ANN has the best performance among the
selected models. In particular, ensemble ANN has the lowest average errors,
$\overline{\text{MAE}}$ = 18 meV and $\overline{\text{RMSE}}$ = 29 meV with
the standard deviation of 1 and 3 meV respectively.

In the $k$-fold CV, the entire dataset is split into $k$ approximately equal
parts (folds). Each ML model is trained on the $k-1$ folds and evaluated on
one fold. The $k$-fold procedure was implemented using the scikit-learn
library scikit-learn  with $k=6$. The results of the $k$-fold CV are presented
in Table 2. Similarly to the shuffle-split CV, $k$-fold CV shows that the
ensemble ANN has the best scores with $\overline{\text{MAE}}$ = 18 meV and
$\overline{\text{RMSE}}$ = 28 meV with the standard deviation of 1 and 3 meV
respectively.

We also evaluated the ANN model on the random test-train split with 20 % of
the data allocated for the test subset. The prediction of transfer integrals
for the test set is shown in Fig. 3 as a scatter plot of the calculated values
versus predicted ones.

##  5 Discussion

In this section, we will apply our ensemble ANN model to different classes of
cuprates and discuss its predictive power. The first example are parent
compounds of high-temperature superconductors (HTSC). The common structural
feature of these antiferromagnets are cuprate planes formed by corner-sharing
CuO4 plaquettes. The Cu-O-Cu angle amounts to 180∘, maximizing electron
transfer between the nearest neighbors ($t_{1}$) and giving rise to a sizable
antiferromagnetic exchange of about 1500 K coldea01 . In addition, the
favorable mutual orientation of plaquettes boosts the coupling between second
neighbors ($t_{2}$), as confirmed experimentally tanaka04 . Hence, the
magnetic properties of undoped HTSC cuprates are described by the frustrated
square-lattice model, with competing first- and second-neighbor
antiferromagnetic exchanges. Interestingly, ramifications of this competition
go far beyond the magnetism: the $t_{2}/t_{1}$ ratio shows correlations with
the superconducting transition temperature Pavarini2001 . Thus, an accurate
estimation of this ratio is crucial for understanding the physics of HTSC
materials.

To test the accuracy of our ensemble ANN model, we consider the following
parent HTSC compounds: La2CuO4, Bi2Sr2CuO6, Tl2Ba2CuO6, HgBa2CuO4,
Tl2Ba2CaCu2O8, HgBa2CaCu2O6, and HgBa2Ca2Cu3O8. It is important to note that
only the former four structures were included in the training dataset. For all
compounds, we recover the frustrated square lattice model with a dominant
$t_{1}$ and a considerably smaller $t_{2}$. For ease of comparison with Ref.
Pavarini2001 , we plot the resulting $t_{2}/t_{1}$ ratios as a function of the
distance between Cu and the apical oxygen atom ($d_{\text{a}}$) in Fig. 4(a).
In the same plot, we show the results of direct calculations of these transfer
integrals by DFT calculations and Wannierization. A very good agreement is
found for all seven cases.

Figure 4: (Color online) (a) The $t_{2}/t_{1}$ ratio as a function of Cu -
apical oxygen distance $d_{a}$ for HTSC cuprates. Black circles corresponds to
the predictions of ensemble ANN model, and orange triangles – results of the
DFT calculations. Numbers correspond to the following compounds:
1-Tl2Ba2Ca2Cu3O10, 2-La2CuO4, 3-Bi2Sr2CuO6, 4-Tl2Ba2CuO6, 5-HgBa2CuO4,
6-Tl2Ba2CaCu2O8, 7-HgBa2CaCu2O6, 8-HgBa2Ca2Cu3O8. (b) Schematics of the
distinct structural element of the HTSC cuprate hosting leading $t_{1}$ and
considerably smaller $t_{2}$ transfer integrals.

The closely related family of double-perovskite cuprates A2CuTO6 (A = Ba or
Sr, T = Te or W) represents a more challenging test case. Here, the
$t_{2}/t_{1}$ ratio crucially depends on the nature of the T atom: in the two
Te-containing compounds, the leading coupling follows the shortest connections
($t_{1}$), while in the other two compounds the empty $5d$ shell of W boosts
the diagonal coupling ($t_{2}$) katukuri20 . The sensitivity of our ANN model
does not suffice to fully account for this trend: it yields a dominant $t_{2}$
for all four compounds. Despite this shortcoming, the model correctly
reproduces the $t_{1}$-$t_{2}$ model, and the predicted $t_{2}/t_{1}$ ratio is
lower for Te-containing (1.55 for Sr2CuTeO6 and 1.3 for Ba2CuTeO6) than for
W-containing (1.6 for Sr2CuTeO6 and 1.95 for Ba2CuWO6) compounds.

Next, we turn to quasi-one-dimensional cuprates. The dominance of $t_{1}$ is
correctly reproduced for the quasi-one-dimensional Sr2CuO3 rosner97 , another
compound with corner-sharing connections between CuO4 squares. Importantly,
this structure features one shorter Cu-Cu connection, which is not accompanied
by sizable electron transfer, and our model correctly captures this aspect:
the respective predicted transfer integral is about 20 times smaller than the
leading intra-chain term. For another quasi-one-dimensional compound, linarite
PbCuSO4(OH)2 featuring edge-sharing chains, our model correctly recognizes the
relevance of first- and second-neighbor transfer integrals along the chains,
and correctly identifies the leading interchain coupling heinze22 .
Importantly, in linarite like in many other edge-sharing cuprates, the
nearest-neighbor exchange is ferromagnetic. Such exchanges have a more complex
nature and can not be described within the effective one-orbital model, which
is at the core of our approach. However, the presence of an edge-sharing
connection essentially implies the relevance of the respective magnetic
exchange, making a dedicated estimation of the electron transfer unnecessary.

As a less trivial case, we consider two isostructural natural minerals in
which Cu2+ atoms form a kagome lattice: kapellasite $\alpha$-Cu3Zn(OH)6Cl2 and
haydeeite $\alpha$-Cu3Mg(OH)6Cl2. The relevance of the cross-hexagon coupling
$t_{d}$ and the corresponding magnetic exchange was suggested based on DFT
results janson08 ; jeschke13  and confirmed experimentally fak12 ; boldrin15 .
(As a side note, the $J_{d}$ exchange is the principal source of frustration
in these systems, because the nearest-neighbor exchange $J_{1}$ is
ferromagnetic.) Here, we consider crystal structures of kapellasite and
haydeeite that were determined by neutron diffraction; these structures are
not in the ICSD and hence were not used for training. The ANN model yields
nearly identical results for both materials, suggesting the leading
$t_{1}\simeq 78$ meV, plus sizable $t_{2}$ and $t_{d}$ of about 25 meV each.
The $t_{1}$ and $t_{d}$ values are comparable with the first-principles
calculations janson08 . Given that the latter used a different functional and
a different structural input, the agreement is very good. Yet, the relevance
of $t_{2}$ in the ensemble ANN model is a spurious result which is at odds
with the DFT calculations and experiments.

In all previous examples except Sr2CuO3, the electronic structure featured
several relevant transfer integrals. There are many cuprates whose magnetism
is shaped by a single coupling dominating over other terms, but it is unclear
which coupling is dominant. An instructive example is the spin-dimer compound
Cu2TeO5, where magnetic dimers do not coincide with the structural ones das08
; ushakov09 . For this material, our ANN model successfully reproduces the
magnitude of the strongest coupling and its dominance over other terms.
Another relevant example is the spin-chain compound CuSe2O5 janson2009 , where
electron transfer is facilitated by the $[$Se2O${}_{5}]^{2-}$ anionic group
connecting two CuO4 squares that are at an angle to each other. Also here the
ANN model correctly identifies the leading transfer integral.

Naturally, the predictive power of the model is limited, and in some cases the
desired accuracy is not reached. For instance, in Bi2CuO4 the leading coupling
operates between the structural chains formed by stacks of CuO4 squares, while
the nearest-neighbor coupling within these stacks is three times weaker
janson07 . Our ANN model correctly reproduces the leading coupling, yet
predicts that the nearest-neighbor coupling has a similar strength. While the
structure of ANN does not allow us to unequivocally determine the root cause
of this discrepancy, we believe that it stems from the correlation between the
magnitude of the transfer integral and the Cu..Cu separation. While on the
average shorter distances indeed correspond to larger $\|t\|$, in some
materials like Bi2CuO4 it is not the case.

There are several ways to improve the accuracy of the model. An apparent
solution is to extend the dataset by including structures that are not
represented in the ICSD. Also a revision of materials that were filtered out
due to failed Wannierization can make the dataset bigger. However, such
amendments will lead to incremental, moderate improvement of the predictive
power. Based on our analysis, we conclude that main factor limiting the
accuracy of the model are the crystalline environments descriptors. Making
them more specific to chemical environments, e.g. by taking the connectivity
of atoms within a chosen sphere into account or a more explicit consideration
of charge densities, and keeping them as compact as possible may significantly
improve the performance of the model.

To finalize the discussion, we emphasize that the main strength our model is
its ability to identify magnetically relevant couplings. This is particularly
useful for involved structures with a large number of short- and middle-range
Cu..Cu separations, where the leading electron transfer paths can be highly
nontrivial. The performance remains good across different classes of cuprates,
which allows for efficient screening: evaluation of transfer integrals for a
single material takes between dozens of seconds and a few minutes. Naturally,
our model can not serve as a complete replacement to full-blown first-
principle calculations, because error bars for the individual terms may be too
high for certain quantitative analyses. However, the model’s predictive power
is enough to perform qualitative assessment of interactions in spin models.
Furthermore, the developed model holds substantial promise for enabling the
inverse construction of hypothetical materials with prescribed magnetic
topologies. For instance, one can create a Cu-O+X network and manipulate its
structure to achieve a particular magnetic coupling arrangement, as indicated
by the transfer integrals data generated by the model. Alternatively, one can
begin with a known material and inquire about the alterations needed to
activate or deactivate, as well as strengthen or weaken, specific magnetic
connections. This then opens a plethora of questions on how these enhanced
properties may be received in an actual chemical solid-state structure.
Thereby this method may offer a novel avenue to engineer materials with
distinct magnetic properties and unlock these applications, for instance in
the fields of magnetic cooling or data storage.

##  6 Conclusions

We constructed an ensemble deep learning model that estimates the magnitude of
transfer integrals in undoped cuprates. These terms underlie the leading
mechanism of the magnetic exchange, and their knowledge is crucial for
correctly determining the microscopic magnetic model. We employed a mapping
onto a three-dimensional Zernike descriptor to describe crystalline
environments that correspond to individual transfer integrals. The resulting
ANN model trained on our high-throughput DFT calculations results can predict
transfer integrals with reasonable error MAE = 18 meV. The model efficiently
differentiates between weak and sizable transfer integrals, which is most
important for estimating the relevant spin model. We discuss the limitations
of this approach and outline ways of improving the numerical accuracy.

Acknowledgments We thank Markus Wallerberger, Roman Rezaev, and Dmitry
Chernyavsky for fruitful discussions. This project was funded by the Leibniz
Association through the Leibniz Competition. Authors acknowledge financial
support from the DFG through the Collaborative Research Center SFB 1143
(Project-Id247310070). We thank Ulrike Nitzsche for the technical assistance.

## Methods

DFT calculations For high-throughput DFT calculations, we used the crystal
structures of cuprates downloaded from the ICSD
bergerhoff.crystallographic.1987  employing the application programming
interface (API). We performed nonmagnetic DFT calculations with the full
potential local orbital code FPLO of version 18.00-52 koepernik.prb.1999 .
Electron exchange-correlation interactions were described using the Perdew-
Burke-Ernzerof (PBE) GGA functional perdew.prl.1996 . The electron density was
converged to within 10-6. The reciprocal space mesh was calculated for each
structure, accounting for the size of the reciprocal cell suppl .

One shot calculation of Hellmann–Feynman forces was performed for each cuprate
structure. We performed the relaxation of structures containing hydrogen by
optimizing atomic coordinates of all H atoms with respect to the GGA total
energy. All symmetries of the respective space group were kept during the
optimization. The Wannier fit of the band structure was performed only for
those structures where calculated forces are below 0.1 eV/Å suppl .

The Python library pymatgen Ong2013pymatgen  has been used broadly in the
high-throughput pipeline code to achieve complete automation.

Machine learning model In the work we consider regression problem
$t_{ij}=f(\vec{D})$ which we handle with ML methods. For selection of the
predictive model we consider shuffle-split and $k$-fold CV. The CV procedures
were implemented with scikit-learn Python library scikit-learn . The best
performance was shown by the ensemble ANN model. We implemented ANN using the
Keras API chollet2015keras  written in Python for learning platform Tensorflow
tensorflow2015-whitepaper . The Supplementary Information 1 suppl  provides
the details on the architecture and training process of the ensemble ANN.

## References

  *   * (1) M. Tsubaki, T. Mizoguchi, Quantum Deep Field: Data-Driven Wave Function, Electron Density Generation, and Atomization Energy Prediction and Extrapolation with Machine Learning.  Phys. Rev. Lett. 125, 206401 (2020).
  * (2) S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine Learning for Fluid Mechanics.  Annu. Rev. Fluid Mech. 52(1), 477–508 (2020).
  * (3) J. Schmidt, M. R. G. Marques, S. Botti, M. A. L. Marques, Recent advances and applications of machine learning in solid-state materials science.  npj Comput. Mater. 5(1), 83 (2019).
  * (4) T. Xie, J. C. Grossman, Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties.  Phys. Rev. Lett. 120, 145301 (2018).
  * (5) M. Rupp, A. Tkatchenko, K. R. Müller, O. A. von Lilienfeld, Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning.  Phys. Rev. Lett. 108, 058301 (2012).
  * (6) K. T. Schütt, H. Glawe, F. Brockherde, A. Sanna, K. R. Müller, E. K. U. Gross, How to represent crystal structures for machine learning: Towards fast prediction of electronic properties.  Phys. Rev. B 89(20), 205118 (2014).
  * (7) A. P. Bartók, R. Kondor, G. Csányi, On representing chemical environments.  Phys. Rev. B 87(18), 184115 (2013).
  * (8) A. Ziletti, D. Kumar, M. Scheffler, L. M. Ghiringhelli, Insightful classification of crystal structures using deep learning.  Nat. Commun. 9(1), 2775 (2018).
  * (9) V. Venkatraman, P. R. Chakravarthy, D. Kihara, Application of 3D Zernike descriptors to shape-based ligand similarity searching.  J. Cheminformatics 1(1), 19 (2009).
  * (10) M. Novotni, R. Klein, in _Proceedings of the eighth ACM symposium on Solid modeling and applications - SM '03_ (ACM Press, 2003).
  * (11) M. Novotni, R. Klein, Shape retrieval using 3D Zernike descriptors.  Comput. Aided Des. 36(11), 1047–1062 (2004).
  * (12) L. Sael, B. Li, D. La, Y. Fang, K. Ramani, R. Rustamov, D. Kihara, Fast protein tertiary structure retrieval based on global surface shape similarity.  Proteins: Structure, Function, and Bioinformatics 72(4), 1259–1273 (2008).
  * (13) L. Mak, S. Grandison, R. J. Morris, An extension of spherical harmonics to region-based rotationally invariant descriptors for molecular shape description and comparison.  J. Mol. Graph. Model. 26(7), 1035–1045 (2008).
  * (14) L. Sael, D. La, B. Li, R. Rustamov, D. Kihara, Rapid comparison of properties on protein surface.  Proteins: Struct. Funct. Genet. 73(1), 1–10 (2008).
  * (15) N. Plakida, _High-temperature cuprate superconductors_ , 2010th edn.  Springer Series in Solid-State Sciences (Springer, Berlin, Germany, 2010).
  * (16) T. Saha-Dasgupta, The Fascinating World of Low-Dimensional Quantum Spin Systems: Ab Initio Modeling.  Molecules 26(6), 1522 (2021).
  * (17) J. S. Helton, K. Matan, M. P. Shores, E. A. Nytko, B. M. Bartlett, Y. Yoshida, Y. Takano, A. Suslov, Y. Qiu, J. H. Chung, D. G. Nocera, Y. S. Lee, Spin Dynamics of the Spin-1/2 Kagome Lattice Antiferromagnet ZnCu3(OH)6Cl2.  Phys. Rev. Lett. 98, 107204 (2007).
  * (18) M. Jaime, V. F. Correa, N. Harrison, C. D. Batista, N. Kawashima, Y. Kazuma, G. A. Jorge, R. Stern, I. Heinmaa, S. A. Zvyagin, Y. Sasago, K. Uchinokura, Magnetic-Field-Induced Condensation of Triplons in Han Purple Pigment BaCuSi2O6.  Phys. Rev. Lett. 93, 087203 (2004).
  * (19) Y. Kohama, H. Ishikawa, A. Matsuo, K. Kindo, N. Shannon, Z. Hiroi, Possible observation of quantum spin-nematic phase in a frustrated magnet.  Proc. Nat. Acad. Sci. 116, 10686–10690 (2019).
  * (20) J. B. Goodenough, _Magnetism And The Chemical Bond_ (John Wiley And Sons, 1963).
  * (21) J. Hubbard, Electron correlations in narrow energy bands.  Proc. R. Soc. Lond. A 276(1365), 238–257 (1963).
  * (22) P. W. Anderson, New Approach to the Theory of Superexchange Interactions.  Phys. Rev. 115, 2–13 (1959).
  * (23) A. Auerbach, F. Berruto, L. Capriotti, _Quantum Magnetism Approaches to Strongly Correlated Electrons_ (Springer Berlin Heidelberg, 2000), pp. 143–170.
  * (24) A. Belik, M. Azuma, M. Takano, Characterization of quasi-one-dimensional S=1/2 Heisenberg antiferromagnets Sr2Cu(PO4)2 and Ba2Cu(PO4)2 with magnetic susceptibility, specific heat, and thermal analysis.  J. Solid State Chem. 177(3), 883–888 (2004).
  * (25) M. D. Johannes, J. Richter, S. L. Drechsler, H. Rosner, ${\mathrm{Sr}}_{2}\mathrm{Cu}{(\mathrm{P}{\mathrm{O}}_{4})}_{2}$: A real material realization of the one-dimensional nearest neighbor Heisenberg chain.  Phys. Rev. B 74, 174435 (2006).
  * (26) O. Janson, W. Schnelle, M. Schmidt, Y. Prots, S. L. Drechsler, S. K. Filatov, H. Rosner, Electronic structure and magnetic properties of the spin-1/2 Heisenberg system CuSe2O5.  New J. Phys. 11(11), 113034 (2009).
  * (27) G. Bergerhoff, I. Brown, F. Allen, et al., Crystallographic databases.  International Union of Crystallography, Chester 360, 77–95 (1987).
  * (28) See Supplementary Information 1 for details of Wannierization, the reciprocal space mesh calculation, and the numerical implementation of functions $\mathcal{I}(x,y,z)$ and $\mathcal{Z}_{nl}^{m}(x,y,z)$. See Supplementary Information 2 for the summary of Wannier fits.
  * (29) S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis.  Comput. Mater. Sci. 68, 314–319 (2013).
  * (30) J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple.  Phys. Rev. Lett. 77(18), 3865–3868 (1996).
  * (31) K. Koepernik, H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme.  Phys. Rev. B 59(3), 1743–1757 (1999).
  * (32) K. Koepernik, O. Janson, Y. Sun, J. van den Brink.  Symmetry Conserving Maximally Projected Wannier Functions (2021).  arXiv:2111.09652.
  * (33) H. Eschrig, K. Koepernik, Tight-binding models for the iron-based superconductors.  Phys. Rev. B 80, 104503 (2009).
  * (34) N. Canterakis, Complete moment invariants and pose determination for orthogonal transformations of 3D objects.  Internal Report 1/96, Technische Informatik I, Technische Universität Hamburg-Harburg (1996).
  * (35) N. Canterakis, 3D Zernike moments and Zernike affine invariants for 3D image analysis and recognition.  11th Scandinavian Conference on Image Analysis (1999).
  * (36) F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, _NIST Handbook of Mathematical Functions_ (Cambridge University Press, 2010).
  * (37) J. Morais, I. Cação, Quaternion Zernike spherical polynomials.  Math. Comput. 84(293), 1317–1337 (2014).
  * (38) L. Breiman, Random forests.  Mach. Learn. 45(1), 5–32 (2001).
  * (39) L. Breiman, Bagging predictors.  Machine Learning 24(2), 123–140 (1996).
  * (40) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine Learning in Python.  J. Mach. Learn. Res. 12, 2825–2830 (2011).
  * (41) R. Coldea, S. M. Hayden, G. Aeppli, T. G. Perring, C. D. Frost, T. E. Mason, S. W. Cheong, Z. Fisk, Spin waves and electronic interactions in ${\mathrm{la}}_{2}{\mathrm{cuo}}_{4}$.  Phys. Rev. Lett. 86, 5377–5380 (2001).
  * (42) K. Tanaka, T. Yoshida, A. Fujimori, D. H. Lu, Z. X. Shen, X. J. Zhou, H. Eisaki, Z. Hussain, S. Uchida, Y. Aiura, K. Ono, T. Sugaya, T. Mizuno, I. Terasaki, Effects of next-nearest-neighbor hopping ${t}^{{}^{\prime}}$ on the electronic structure of cuprate superconductors.  Phys. Rev. B 70, 092503 (2004).
  * (43) E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, O. K. Andersen, Band-structure trend in hole-doped cuprates and correlation with ${\mathit{t}}_{\mathit{c}\mathrm{max}}$.  Physical Review Letters 87(4) (2001).
  * (44) V. M. Katukuri, P. Babkevich, O. Mustonen, H. C. Walker, B. Fåk, S. Vasala, M. Karppinen, H. M. Rønnow, O. V. Yazyev, Exchange interactions mediated by nonmagnetic cations in double perovskites.  Phys. Rev. Lett. 124, 077202 (2020).
  * (45) H. Rosner, H. Eschrig, R. Hayn, S. L. Drechsler, J. Málek, Electronic structure and magnetic properties of the linear chain cuprates ${\mathrm{Sr}}_{2}{\mathrm{CuO}}_{3}\phantom{\rule{0.0pt}{0.0pt}}$ and ${\mathrm{Ca}}_{2}{\mathrm{CuO}}_{3}$.  Phys. Rev. B 56, 3402–3412 (1997).
  * (46) L. Heinze, M. D. Le, O. Janson, S. Nishimoto, A. U. B. Wolter, S. Süllow, K. C. Rule, Low-energy spin excitations of the frustrated ferromagnetic $J_{1}-J_{2}$ chain material linarite PbCuSO4(OH)2 in applied magnetic fields parallel to the $b$ axis.  Phys. Rev. B 106, 144409 (2022).
  * (47) O. Janson, J. Richter, H. Rosner, Modified Kagome Physics in the Natural Spin-1/2 Kagome Lattice Systems: kapellasite Cu3Zn(OH)6Cl2 and Haydeeite Cu3Mg(OH)6Cl2.  Phys. Rev. Lett. 101, 106403 (2008).
  * (48) H. O. Jeschke, F. Salvat-Pujol, R. Valentí, First-principles determination of Heisenberg Hamiltonian parameters for the spin-$\frac{1}{2}$ kagome antiferromagnet ZnCu3(OH)6Cl2.  Phys. Rev. B 88, 075106 (2013).
  * (49) B. Fåk, E. Kermarrec, L. Messio, B. Bernu, C. Lhuillier, F. Bert, P. Mendels, B. Koteswararao, F. Bouquet, J. Ollivier, A. D. Hillier, A. Amato, R. H. Colman, A. S. Wills, Kapellasite: A Kagome Quantum Spin Liquid with Competing Interactions.  Phys. Rev. Lett. 109, 037208 (2012).
  * (50) D. Boldrin, B. Fåk, M. Enderle, S. Bieri, J. Ollivier, S. Rols, P. Manuel, A. S. Wills, Haydeeite: A spin-$\frac{1}{2}$ kagome ferromagnet.  Phys. Rev. B 91, 220408 (2015).
  * (51) H. Das, T. Saha-Dasgupta, C. Gros, R. Valentí, Proposed low-energy model Hamiltonian for the spin-gapped system ${\text{CuTe}}_{2}{\text{O}}_{5}$.  Phys. Rev. B 77, 224437 (2008).
  * (52) A. V. Ushakov, S. V. Streltsov, Electronic and magnetic structure for the spin-gapped system CuTe2O5.  J. Phys.: Condens. Matter 21, 305501 (2009).
  * (53) O. Janson, R. O. Kuzian, S. L. Drechsler, H. Rosner, Electronic structure and magnetic properties of the spin-1/2 Heisenberg magnet Bi2CuO4.  Phys. Rev. B 76, 115119 (2007).
  * (54) F. Chollet, et al.  Keras.  https://keras.io (2015).
  * (55) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems (2015).  Software available from tensorflow.org.


11institutetext: Univ. Lille, ULR 7512 - Unité de Mécanique de Lille Joseph
Boussinesq (UML), F-59000 Lille, France

# Surface light modulation by sea ice and phytoplankton survival in a
convective flow model

Vinicius Beltram Tergolina  Enrico Calzavarini  Gilmar Mompean  Stefano Berti

(Received: date / Revised version: date)

###### Abstract

Plankton dynamics depend in a complex manner on a variety of physical
phenomena, according to both experimental and numerical data. In particular,
experimental field studies have highlighted the relation between phytoplankton
survival and turbulent upwelling and downwelling from thermal convection.
Recent numerical works have also shown the importance of accounting for
advective transport by persistent structures in simulation models. In
nutrient-rich polar marine environments phytoplankton blooms are critically
limited by light availability under ice-covered waters. Such heterogeneity of
the light intensity distribution, in association with a large-scale coherent
fluid flow, can give rise to nontrivial growth dynamics. In this work we
extend a previous advection-reaction-diffusion model of phytoplankton light-
limited vertical dynamics in the presence of convective transport.
Specifically, we consider horizontally heterogeneous light conditions through
the use of two regions with different production regimes, modelling the
absence (presence) of light under (in between) obstacles. Such a model is
intended as an idealized representation of nonuniformly ice-covered polar
waters. By means of numerical simulations, we find that the main role of
advective transport is to hinder phytoplankton growth, but also that such
effect depends on the positions of the obstacles with respect to the upwelling
and downwelling flow regions. Furthermore, we show that the sinking speed due
to the density difference between phytoplankton organisms and water, while
small, plays an important role, which depends on how it adds to the flow.
These results indicate that advective transport can have a crucial impact on
the survival conditions of sinking phytoplankton species in polar
environments.

###### pacs:

PACS-keydiscribing text of that key and PACS-keydiscribing text of that key

##  1 Introduction

In phytoplankton ecology, biological and fluid dynamics are closely linked,
ruling how planktonic populations interact with their environment and other
organisms. Interdisciplinary work aimed at studying the interplay between
plankton biology and the concurrent physical factors has allowed to
considerably refine the understanding of the subject ML2005 ;
kiorboe2018mechanistic ; prairie2012biophysical . Notably, both modeling
studies (see, e.g., taylor2011shutdown ; lindemann2017dynamics ;
tergolina2021effects ) and direct measurements lowry2018under  demonstrated
the importance of vertical advective transport, underlining the relation
between turbulent convective fluid motions and phytoplankton growth dynamics.

Polar marine ecosystems are the seat of important plankton blooms (see
ardyna2020under  for a review). Understanding the dynamics of the latter is
key to correctly estimate primary productivity lowry2018under ;
arrigo2012massive ; arrigo2015continued ; assmy2017leads , particularly under
a changing climate scenario. In the Arctic ocean, several studies reported
under-ice phytoplankton blooms that were initiated by the onset of ice melt
lowry2018under ; arrigo2012massive ; assmy2017leads ; boles2020under .
Nevertheless, it is still not known whether under-ice blooms are controlled by
the increase of light transmittance, leads (openings in the ice), upwelling
from springtime convective mixing, biological factors, or other factors (most
likely a combination of more than one). Interestingly, however, it was
observed that waters beneath sea ice with leads had significantly lower
phytoplankton biomass, despite high nutrient availability, than regions
without (or with few) leads lowry2018under . It was argued that this
phenomenon is controlled by the onset of convection below leads, which
enhances mixing to larger depths, and reduces the travel time along the
vertical. Consistently with Sverdrup’s arguments sverdrup1953conditions ,
then, the bloom should be suppressed due to the mixed-layer depth being larger
than the critical depth (where the integrated production is balanced by the
integrated loss). Note that convection is here to be intended both as a
thermo-solutal buoyancy-induced phenomenon, due to brine rejection from the
ice and surface heat losses, and as wind-induced mixing in leads
lowry2018under . It is worth remarking, too, that phytoplankton blooms beneath
loosely consolidated ice will likely become more common in the future Arctic
ocean, as lead formation becomes more frequent due to increasingly thinner and
dynamic ice coverage lindsay2005thinning ; serreze2007perspectives ;
kwok2013arctic . This could significantly alter primary production, and have
important consequences on local marine food webs.

In this work we take an approach based on idealized models, to explore the
effect on phytoplankton dynamics of the horizontal modulation of light
intensity at the surface of a vertical fluid layer. For this purpose, we
extend a previous two-dimensional (2D) advection-reaction-diffusion model that
we used to investigate the role of large and small scale fluid motions on the
same biological system, but with uniform surface light intensity
tergolina2021effects . The role of vertical turbulent diffusivity on light-
limited phytoplankton dynamics has been extensively studied in the past
shigesada1981analysis ; huisman1994light ; huisman2002sinking ;
huisman2002population  in one-dimensional (1D) systems. Our model can be seen
as a 2D generalization of the one introduced in huisman2002sinking . Moreover,
we also consider the presence of a coherent large-scale flow (as in
tergolina2021effects ) and of horizontally varying incident light intensity at
the surface. The flow field will be given in terms of a streamfunction
corresponding to a kinematic cellular flow. While this is clearly a strong
simplification with respect to a truly dynamic velocity field, it allows us to
use realistic biological parameters and domain size. The cellular flow
structure is intended as a simplified description of buoyancy and wind-driven
flows thorpe2005turbulent , such as convective motions and Langmuir
circulations craik1976rational ; bees1997planktonic ; thorpe2004langmuir ,
often encountered in the upper ocean. Beyond the environmental implications
discussed above, from a more fundamental point of view, based on the
overwhelming importance of the large-scale flow and of its downwelling region
that we found in a previous work tergolina2021effects , in this study we
wanted to see how phyotplankton dynamics change when reproduction is forbidden
over either the upwelling or the downwelling flow region. As also in this case
we expect a weak effect of small-scale flow features, we only consider a
steady large-scale (cellular) flow. We then place a light-blocking obstacle,
either over the downwelling region or over the upwelling one.

Through this approach we aim at assessing the importance of the different
transport mechanisms on phytoplankton survival, in the presence of such
vertically and horizontally heterogeneous growth dynamics. We also examine
different values of the sinking speed, due to phytoplankton being typically
heavier than water. This corresponds to considering different algal species,
and appears interesting on the basis of the experimental evidences about the
varied species composition of phytoplankton in under-ice blooms
arrigo2012massive .

This article is organized as follows. The model is introduced in Sec. 2. In
Sec. 3 we report the results of our numerical study, and we discuss the impact
of the different transport terms, by separately considering them in a
systematic way. The main outcome and the conclusions are presented in Sec. 4.

##  2 Model

Our model is intended as an idealized representation of phytoplankton dynamics
in waters that are partially covered by ice. Specifically, we consider light-
limited growth in a 2D vertical fluid layer, in the presence of obstacles at
the surface mimicking ice blocks that do not allow the transmission of the
incoming radiation to the water underneath. A schematic is provided in Fig. 1.
While our description is an over-simplification of a realistic under-ice
ecosystem, the present setup can provide indications useful to develop more
complex biophysical models, and is in principle quite easily reproducible in
laboratory experiments.

To describe the dynamics of the system is then necessary to specify the
evolution equation of the population density field $\theta(x,z,t)$ (number of
individuals per unit volume) in both the ice-free regions and in the ice-
covered ones. In the regions receiving light (where phytoplankton can
reproduce), we have:

$\frac{\partial\theta}{\partial
t}=\left[p(I)-l\right]\theta-\bm{v}\cdot\bm{\nabla}\theta+D\bm{\nabla}^{2}\theta,$
(1)

as in tergolina2021effects , where $p(I)$ is the light-dependent production
rate and $l$ the loss rate. In the regions below obstacles, instead, due to
the absence of light, $p(I)=0$ and one has:

$\frac{\partial\theta}{\partial
t}=-l\theta-\bm{v}\cdot\bm{\nabla}\theta+D\bm{\nabla}^{2}\theta.$ (2)

The fluid layer in which we examine these dynamics is intended to represent
the mixed layer (the water layer of nearly constant turbulent intensity found
in upper oceans and lakes vallis2017atmospheric ). Its horizontal and vertical
sizes will be denoted $L_{x}$ and $L_{z}$, respectively. The top and bottom
boundaries are assumed to be rigid (as in tergolina2021effects ;
huisman2002sinking ).

Figure 1: Schematic of the model setup. The vertical decrease of light
intensity from the surface value $I_{in}$ is described by Eq.
(4).Phytoplankton (green dots) dynamics are assumed to evolve according to
Eqs. (1) [with (3) and (4)] and (2), in the water columns in between the
surface obstacles and below them (where $I_{in}=0$), respectively. A cellular
flow characterized by both upwelling and downwelling regions is also present.

The production term in Eq. (1) accounts for both water background turbidity,
via the coefficient $\kappa_{bg}$, and population self-shading, with an
attenuation factor $\kappa$. Its functional form is:

$p(I)=\frac{p_{max}I}{H+I},$ (3)

where $p_{max}$ is the maximum specific production rate, $H$ is a half-
saturation constant and the time- and depth-dependent light-intensity is
expressed as follows, according to Lambert-Beer’s law:

$I(z,t)=I_{in}e^{-\int_{0}^{z}\kappa\theta(s,t)ds-\kappa_{bg}z},$ (4)

with $I_{in}$ the incident light (at the surface, where $z=0$). Note that we
will assume $I_{in}=0$ in correspondence with the surface obstacles. The
biological parameter values used in our study, representative of realistic
situations huisman1999critical , are reported in Table 1. They correspond to
average measured values, with growth parameters obtained from freshwater
phytoplankton species, and $\kappa_{bg}$ from clear lakes and coastal areas
huisman2002population .

Table 1: Parameters of the biological dynamics. Parameter | Value | Meaning
---|---|---
$\kappa_{bg}$ |  $0.2$ m-1 | Background turbidity
$\kappa$ |  $1.5\cdot 10^{-11}$ m2 cell-1 | Specific light attenuation of phytoplankton
$p_{max}$ |  $0.04$ h-1 | Maximal specific production rate
$l$ |  $0.01$ h-1 | Specific loss rate
$H$ |  $30\,\mu$mol photons m-2 s-1 | Half-saturation constant of light-limited growth
$I_{in}$ |  $350\,\mu$mol photons m-2 s-1 | Incident light intensity
$v_{sink}$ |  $0.04$ m h-1 | Phytoplankton sinking velocity

The coefficient $D$ in Eqs. (1) and (2) represents an effective diffusivity,
due to both small-scale unresolved turbulent motions and possible non-directed
swimming behavior. Advection is realized by a 2D incompressible flow
$\bm{u}=(u_{x},u_{z})$ and phytoplankton is assumed to sink with a speed
$v_{sink}\hat{\bm{z}}$, where $\hat{\bm{z}}$ is the unitary vector pointing
downward in the vertical direction; the total velocity appearing in the
evolution equations is thus $\bm{v}=\bm{u}+v_{sink}\hat{\bm{z}}$. For the
phytoplankton sinking speed we take the value in Table 1 as the reference one,
and we will later vary it.

For the velocity field we choose a cellular flow, which was originally
introduced as a kinematic model of thermal convection fluid motions
solomon1988chaotic . Its expression is obtained as
$\bm{u}=\left(-\partial_{z}\Psi,\partial_{x}\Psi\right)$ from the following
streamfunction:

$\Psi(x,z)=-\frac{U}{k}\sin(kx)\sin(kz),$ (5)

where $U$ is the maximum velocity and $k=\pi/L_{z}$ a wavenumber. The flow
pattern then corresponds to two counter-rotating steady vortices of size
$L_{z}$, in the $(x,z)$ plane. It is worth stressing that this flow comprises
both a (central) downwelling region and two (side) upwelling ones, which
typically play a relevant role on phytoplankton dynamics lindemann2017dynamics
; tergolina2021effects . Note that eddying motions on scales smaller than
$L_{z}$, and unsteadiness, could be in principle taken into account by adding
to the streamfunction of Eq. (5) analogous terms with larger wavenumbers, and
an explicit time dependency (see, e.g., tergolina2021effects ;
lacorata2017chaotic ; forgia2022numerical ). Our choice to neglect them, here,
is motivated by a previous study without the surface light modulation
tergolina2021effects , where we found that small unsteady eddies play a
secondary role. Indeed, in that case, the fate (survival or extinction) of
phytoplankton was primarily determined by the large-scale coherent component
of the flow.

In the following we will consider two different geometrical configurations,
corresponding to light-blocking obstacles, of fixed total size, positioned
either above the downwelling flow region (LBD case), as in Fig. 1, or above
the upwelling one (LBU case). The horizontal and vertical domain sizes will be
respectively set to $L_{x}=60$ m and $L_{z}=30$ m. Considering the chosen
biological parameters (Table 1), a mixed layer depth of $30$ m is smaller than
the critical depth (about $50$ m). Hence, in the absence of surface light
heterogeneity and of a fluid flow, with the values of diffusivity explored,
$1$ cm2 s${}^{-1}\leq D\leq 20$ cm2 s-1, which appear quite realistic for the
ocean tergolina2021effects ; huisman2002sinking ; huisman1999critical , a
bloom occurs. The presence of obstacles considerably reduces production, due
to the absence of light below them. However, we made sure that also in this
case phytoplankton blooms when $U=0$. This allowed us to investigate the
impact of advection on the dynamics of the system, including a possible
transition from population survival to extinction. Remark, too, that working
with a relatively small system, but still apt to capture the salient dynamics,
is computationally convenient.

We numerically study the dynamics of Eqs. (1) and (2), which are
simultaneously integrated in our rectangular domain (with $L_{x}=2L_{z}$) by
means of a pseudo-Lagrangian algorithm abel2001front ; berti2005mixing  (see
tergolina2021effects  for more details). We adopt periodic and no-flux
boundary conditions along the horizontal ($x$) and the vertical ($z$),
respectively (as in tergolina2021effects ). The initial condition is a low
uniform population density [$\theta(t=0)=5.5\cdot 10^{6}$ cells m-3]. To
analyze the blooming conditions of the phytoplankton population we rely on the
temporal behavior of the average biomass density,

$\langle\theta\rangle(t)=\frac{1}{L_{x}L_{z}}\int_{0}^{L_{x}}\int_{0}^{L_{z}}\theta(x,z,t)\,dxdz,$
(6)

and the per-capita growth rate,

$r_{p}(t)=\frac{1}{\langle\theta\rangle}\frac{\partial\langle\theta\rangle}{\partial
t}.$ (7)

The per-capita growth rate is expected to attain a statistically constant
value $r_{p}$ (after an initial transient), corresponding to exponential
growth ($r_{p}>0$) or decay ($r_{p}<0$) in the early regime before the onset
of nonlinear dynamical effects.

##  3 Results

Figure 2: Instantaneous normalized population density field
$\theta(x,z,t^{*})/\langle\theta\rangle$ at a fixed instant of time
$t^{*}=600$ h, for $D=10$ cm2 s-1 and $U=(0,0.71,1.19)$ m h-1 [panels (a,d),
(b,e) and (c,f), respectively], where $\langle\theta\rangle$ stands for the
spatial average. The obstacle is represented by the black rectangle, the
vertical dashed lines delimit the light obstruction zone. The white lines are
the isolines $\theta/\langle\theta\rangle=1$, and $t^{*}$ is in the regime of
stationary per-capita growth rate ($r_{p}(t)=r_{p}=\mathrm{const}$). The solid
black lines represent flow streamlines, with arrows indicating the circulation
direction.

We begin our numerical study by exploring the dynamics in two partially
illuminated setups, both in the absence and in the presence of the fluid flow.
As explained in Sec. 2, the blocking of light at the surface is realized by
introducing obstacles of zero light transmittance, with different positions
according to the cases (LBD, LBU) considered. Given the large number of
control parameters (for both biology and transport) in the model, we decided
to work with a single light-transmittance value associated with the obstacle.
We then chose to use fully absorbing obstacles, in order to increase the
contrast between the zones where light availability permits reproduction and
those where this does not occur. Such a simplifying choice is also intended to
ease the interpretation of the effects of the different transport and
production regimes. Indeed, as we shall see, already in the present case the
dynamics are quite complex. The light transmittance of ice is nevertheless
variable and is a subject of debate (see nicolaus2012changes  for recent
measurements of ice thinning in the Arctic ocean).

Some snapshots of the field $\theta$ at a given time are shown in Fig. 2, for
both the LBD and LBU cases, and for increasing values of the flow intensity
$U$. This figure is obtained with the reference sinking speed $v_{sink}=0.04$
m h-1 tergolina2021effects ; huisman2002sinking . When advection is not
present the two light-modulation cases give identical results, as expected.
With nonzero advection, instead, important differences appear. While in the
LBD case phytoplankton accumulates closer to the surface and presents small
differences with respect to increasing flow intensity, in the LBU case the
population spreads more in the vertical, reaching larger depths. By measuring
the per-capita growth rate $r_{p}$ versus time for different values of $U$, we
could assess that phytoplankton growth is hindered by advection, independently
of the light-modulation considered. Indeed, the (constant) value of $r_{p}$
reached at large enough times, monotonically decreases with increasing $U$ in
both the LBD and LBU cases [see Fig. 3, where $r_{p}$ is normalized by the
total (birth minus death) growth rate at the surface
$r_{b}=p_{max}I_{in}/(H+I_{in})-l$]. This effect parallels the one we already
found in the absence of obstacles tergolina2021effects . However, the decrease
of $r_{p}$ in the two situations is quantitatively different, and indeed
larger when light does not penetrate in the upwelling region, so that an
LBD/LBU asymmetry emerges. In the remainder of this work, we focus on the
understanding of this feature.

Figure 3: Population per-capita growth rate $r_{p}(t)$, normalized by the
intrinsic total (surface) growth rate $r_{b}$, as a function of time, for
various values of $U$, $D=10$ cm2 s-1 and $L_{z}=30$ m. Panels (a) and (b)
respectively refer to the LBD and LBU cases.

We can first observe that the orientation, with respect to the light gradient,
of the fluid velocity alone (which is opposite in the LBD and LBU cases)
cannot be responsible for this effect. To see this, one can rewrite Eq. (1) as
follows:

$\frac{\partial\theta}{\partial
t}=g(\bm{x})\theta-\bm{u}\cdot\bm{\nabla}\theta,$ (8)

where $g(\bm{x})$ indicates the spatially dependent growth rate, and diffusion
and sinking have been neglected. Note that in our system the growth rate
depends on both the vertical and the horizontal coordinates, due to the
external light blocking and to self-shading. The fluid velocity
$\bm{u}(\bm{x})$ depends on space, too, but not on time. In a Lagrangian
reference frame the above equation reads:

$\displaystyle\begin{cases}\dot{\theta}=g(\bm{x}(t))\theta,\\\
\dot{\bm{x}}(t)=\bm{u}(\bm{x}(t)).\end{cases}$ (9)

Given an intial value for $\theta$ at time $t=0$ and a final position
$\bm{x}(t)$, the solution of Eq. (9) is:

$\displaystyle\theta(\bm{x}(t),t)$ $\displaystyle=$
$\displaystyle\theta(\bm{x}(0),0)\
e^{\int_{0}^{t}g(\bm{x}(t^{\prime}))dt^{\prime}}$ (10) $\displaystyle=$
$\displaystyle\theta(\bm{x}(0),0)\
e^{\int_{|\bm{x}(0)|}^{|\bm{x}(t)|}\frac{g(\bm{x}^{\prime})}{|\bm{u}(\bm{x}^{\prime})|}dx^{\prime}}.$
(11)

It is therefore clear that the growth of $\theta$ along a given Lagrangian
trajectory (coinciding with a streamline, due to flow stationarity) depends on
the velocity at which the path is travelled but not on the orientation of the
velocity itself. Moreover along a closed path, like a cellular-flow
streamline, the variation of the reactive-scalar density will be independent
of the circulation direction. A full loop starting with downward flow and
followed by an upward flow will produce the same phytoplankton growth if it
occurs in the reversed order.

The presence of a sinking speed, which modifies the Lagrangian equation of
motion to $\dot{\bm{x}}(t)=\bm{u}(\bm{x}(t))+v_{sink}\hat{\bm{z}}$, can change
the scenario because a descending flow path will be covered faster than an
ascending one, leading to a smaller growth for the first than for the latter.
Note that the presence of sinking also has an effect on the shapes of
planktonic trajectories, in particular not all the trajectories in the
cellular flow will be closed stommel1949trajectories . It is therefore the
(slow) sinking movement that is responsible for breaking the symmetry of the
growth dynamics between the upwelling and downwelling flow regions. The
presence of diffusion will also act in a similar manner by allowing parcels of
population density to cross the flow streamlines.

Figure 4 reports the measure of the per-capita growth rate - normalized by the
total (birth minus death) growth rate at the surface $r_{b}$ \- versus time,
for different values of the sinking speed and a fixed flow intensity. Panel
(a) in the figure refers to a small value of the diffusion coefficient, while
panel (b) is for a larger one. The role of diffusivity is quite subtle and we
will come back to it later in this section. Independently of this, the results
clearly show that in the absence of sinking ($v_{sink}=0$) the system is
insensitive to the position of the obstacle. The value of $r_{p}$ is in fact
the same, within numerical accuracy, in the LBD and LBU cases. However, it
decreases in the presence of sinking when light is blocked over the upwelling
region, as it may be expected. It is perhaps more surprising that it instead
increases in the LBD case (which is more difficult to see for large $D$ in
Fig. 4b due to the scale). For small diffusivity, it can even turn a no-bloom
situation ($r_{p}<0$) into a bloom one ($r_{p}>0$).

Figure 4: Population per-capita growth rate, normalized by the intrinsic total
(birth minus death) growth rate at the surface, as a function of time for
different values of the sinking speed in both the LBD case and the LBU one.
The flow intensity is $U=0.936$ m/h, the diffusion coefficient is either $D=1$
cm2 s-1 (a) or $D=10$ cm2 s-1 (b).

The combined effect of fluid advection, sinking and diffusion can be
appreciated through the following argument, which helps rationalizing the
previous observations. We consider the transformation (see also
ryabov2008population ):

$\theta(\bm{x},t)=\textrm{exp}\left(\frac{\bm{w}\cdot\bm{x}}{2D}\right)\tilde{\theta}(\bm{x},t)$
(12)

where $\bm{w}$ is a, yet to be defined, constant (in both time and space)
velocity. Plugging Eq. (12) into Eq. (1) we get:

$\frac{\partial\tilde{\theta}}{\partial
t}=\left[p(I)-l-\frac{\bm{w}\cdot\bm{v}}{2D}+\frac{w^{2}}{4D}\right]\tilde{\theta}-(\bm{v}-\bm{w})\cdot\bm{\nabla}\tilde{\theta}+D\bm{\nabla}^{2}\tilde{\theta}.$
(13)

Taking $\bm{w}=v_{sink}\hat{\bm{z}}$, and using the fact that
$\bm{v}=\bm{u}+v_{sink}\hat{\bm{z}}$, we obtain

$\frac{\partial\tilde{\theta}}{\partial
t}=\left[p(I)-l-\frac{v_{sink}u_{z}}{2D}-\frac{v_{sink}^{2}}{4D}\right]\tilde{\theta}-\bm{u}\cdot\bm{\nabla}\tilde{\theta}+D\bm{\nabla}^{2}\tilde{\theta}.$
(14)

This is the evolution equation of a field advected by the fluid flow $\bm{u}$,
where sinking is not explicitly present, but with a modified reaction term
that now depends on the sign of the vertical velocity component $u_{z}$ (and
also on space). It is thus clear that the local growth rate, i.e.
$g=p(I)-l-v_{sink}u_{z}/(2D)-v_{sink}^{2}/(4D)$, is always smaller in the
downwelling regions ($u_{z}>0$) than in the upwelling ones ($u_{z}<0$). In
short, $g_{d}<g_{u}$, where $d,u$ label down- and upwelling regions,
respectively. We also note that the quadratic term in $v_{sink}$ can be
neglected nearly everywhere, because in practice $v_{sink}\ll|u_{z}|\sim U$
(this may be not true close to the top and bottom boundaries, where $u_{z}$
also becomes very small). The above reasoning is independent of the presence
of the obstacle blocking light at the surface. When the latter is present we
can expect that the LBU situation is characterized by a growth rate similar to
the one observed in downwelling regions, $g_{LBU}\sim g_{d}$ and, on the
opposite, for the LBD case: $g_{LBD}\sim g_{u}$. This implies that
$g_{LBD}>g_{LBU}$, in agreement with the observed trends.

One may then be tempted to extrapolate this result to the per-capita growth
rate, $r_{p}=\partial_{t}\langle\theta\rangle/\langle\theta\rangle$, to obtain
the following predictions for the two situations considered:

$\displaystyle r_{p,LBD}$ $\displaystyle\sim$ $\displaystyle\langle
p(I)\rangle_{L}-l+\frac{2v_{sink}U}{\pi^{2}D}-\frac{v_{sink}^{2}}{4D},$ (15)
$\displaystyle r_{p,LBU}$ $\displaystyle\sim$ $\displaystyle\langle
p(I)\rangle_{L}-l-\frac{2v_{sink}U}{\pi^{2}D}-\frac{v_{sink}^{2}}{4D},$ (16)

which, however, assume a weak spatial coupling between the growth rate $g$ and
$\tilde{\theta}$, and where we used the average vertical velocity in the
regions receiving light, $\langle u_{z}\rangle_{L}=\pm(2/\pi)^{2}U$. More
importantly, here we have neglected any other effects due to both advection
and diffusion. In particular, this argument completely ignores the lateral
transport, between the illuminated and shaded regions, due to both advection
and diffusion. The flux associated with this process is
$u_{x}\theta-D\partial_{x}\theta$, with $u_{x}=\pm Ucos(\pi z/L_{z})$ at the
boundaries between the two regions, and it should be included to improve the
present prediction. Indeed, the expressions in Eqs. (15) and (16) fail to
describe the behavior of the large-time value of $r_{p}$ as a function of
$v_{sink}$, unless some fitting coefficients are introduced to account for the
missing physics. While for small diffusivity the fitting coefficients are
comparable in the LBD and LBU cases, which may suggest that the missing
information resides in the $z$-dependence of $u_{z}$ and in the appropriate
choice of the region where to average, this is not the case for large $D$. In
fact, in this case, a striking feature that is already apparent from Fig. 4b
is the asymmetric response of the system with respect to the position of the
obstacle. This fact cannot be explained by Eqs. (15) and (16), which are
symmetric under the exchange of $U$ into $-U$.

Taking into account the role of lateral transport between shaded and
illuminated regions analytically is not trivial. Remark, too, that in the
large $D$ limit analytical predictions are not available even in the 1D no-
flow (fully illuminated) model huisman2002sinking , starting from which ours
was constructed. In the following we then resort to an argument based on
Lagrangian paths to provide a qualitative explanation of the asymmetric
response of the system.

Figure 5: Distribution of $N=2000$ particles at a time $t>0$ for $D=1$ cm2 s-1
(red circles) and $D=10$ cm2 s-1 (blue squares) in the LBD (a) and LBU (b)
cases. In each case, at time $t=0$ all particles occupied a small area close
to the top of the region receiving light (green circles).

From a Lagrangian point of view, in the presence of diffusion, the complete
equation of motion describing fluid transport is
$\dot{\bm{x}}=\bm{u}(\bm{x}(t))+v_{sink}\hat{\bm{z}}+\sqrt{2D}\bm{\eta}$,
where the white noise term $\sqrt{2D}\bm{\eta}$ represents diffusion. Sample
distributions, at an arbitrary time $t>0$, of $N=2000$ particles evolving
according to this equation are shown in Fig. 5 for $D=1$ cm2 s-1 (red circular
points) and $D=10$ cm2 s-1 (blue square points). The two panels correspond to
the LBD (Fig. 5a) and the LBU (Fig. 5b) case. In each of them, the particles
were initially located in a small area close to the top of the illuminated
region (green points). As expected, for larger $D$ particles spread more into
the domain. It is here worth emphasizing the lateral spreading: particles
originally in a well-lit region are much more likely to reach a dark region,
where reproduction is not possible, if their diffusivity is larger. In the LBD
case this implies that, when $D$ is large enough, lateral diffusion can
overcome the action of the sinking velocity, which helped maintaining
particles in the upwelling reproductive region. In the LBU case, particles are
also more easily brought to deadly regions by increased diffusivity, but this
occurs while sinking contributes moving them downwards, where growth is also
importantly reduced due to light absorption along the vertical. Therefore, the
overall effect is weaker than in the previous case. This explains the
asymmetry observed in the behavior of the per-capita growth rate as a function
of the sinking speed, with respect to the position of the obstacle blocking
light. Indeed, in the LBD case, the positive action of $v_{sink}$ for growth
is only barely visible for $D=10$ cm2 s-1 in Fig. 4b, as it is masked by the
impact of diffusivity. On the contrary, in the LBU case, such a large value of
the diffusion coefficient does not change the general phenomenology (with
respect to the $D=1$ cm2 s-1 case of Fig. 4a), as this was already unfavorable
for growth, due to downward advection by both the sinking and vertical fluid
velocities in the illuminated downwelling region.

We conclude this discussion by noting that the same reasoning can be made also
considering transport backwards in time, which allows interpreting the results
in terms of the generalization of Eq. (10) accounting also for diffusive
transport:

$\theta(\bm{x},t)=\langle\theta(\bm{x}(0),0)\
e^{\int_{0}^{t}g(\bm{x}(t^{\prime}))dt^{\prime}}\rangle_{\eta}\,.$ (17)

The above equation (see abel2001front  for more details) states that the
population density at time $t>0$ and at position $\bm{x}$ is constructed by
advancing in time the reactive dynamics, starting from the value of the field
$\theta(\bm{x}(0),0)$ at the Lagrangian origin $\bm{x}(0)$ of a trajectory
ending at $\bm{x}$ at time $t$, and averaging over noise ($\eta$)
realizations. Clearly, when the diffusion coefficient increases many more
Lagrangian paths reaching $\bm{x}$ at time $t$ originate from regions where
light cannot penetrate, bringing already dead material at the considered
position at a later time. In the LBD case, the effect is considerable, as the
spatially and temporally varying per-capita growth rate $g(\bm{x}(t))$ in Eq.
(17), at a location $\bm{x}$ in the illuminated (and upwelling) region, is now
determined by the many trajectories coming from dark regions, which was not
the case for smaller $D$. In the LBU case, instead, also with small values of
$D$ the contributions to $g(\bm{x}(t))$, with $\bm{x}$ in the illuminated (and
downwelling) region, are essentially associated with trajectories starting in
regions where light is blocked and the population has already died. Therefore,
increasing $D$ does not alter the qualitative behavior of the system in this
case.

##  4 Conclusions

We numerically investigated the dynamics of sinking phytoplankton in a stirred
fluid layer in which light availability is spatially inhomogeneous, along both
the vertical and the horizontal. Our system is assumed to be illuminated from
above, as is typical in aquatic ecosystems receiving energy from the incoming
solar radiation. Vertically, light intensity decreases due to absorption by
water and phytoplankton self-shading. Horizontally, its penetration from the
surface into the interior of the medium is critically affected by the presence
of opaque obstacles at the top boundary.

We focused on the impact of advection, and more generally of the different
transport processes occurring in the fluid, on the average biomass. In
particular, we considered a 2D kinematic velocity field that provides an
idealized representation of the flow rolls observed in thermal convection
solomon1988chaotic , or in Langmuir circulation craik1976rational ;
bees1997planktonic ; thorpe2004langmuir  in the ocean mixed layer. Such flow
pattern is characterized by both upwelling and downwelling regions, which play
a key role on light-limited phytoplankton growth dynamics taylor2011shutdown ;
lindemann2017dynamics ; tergolina2021effects ; lowry2018under . Indeed, a
general feature emerging from previous studies (taking into account the
vertical decrease of light intensity but not its horizontal variability) is
that vertical mixing due to convective motions reduces the growth of
phytoplankton, and may even lead to extinction events lindemann2017dynamics ;
tergolina2021effects . This is largely due to the accumulation of the
population in the downwelling regions associated with persistent large-scale
vortical structures. Consequently, phytoplankton is efficiently transported to
depth, where reproduction is inhibited by the lack of light
lindemann2017dynamics ; tergolina2021effects . Irregular fluid motions on
smaller time and length scales only weakly affect this phenomenology
tergolina2021effects .

Our goal in this work was to explore if and how the response of the system, in
terms of average biomass and survival conditions, changes when the incident
light is not uniform at the surface, a situation that is often encountered in
natural environments. A remarkable instance of such a configuration is found
in oceanic waters that are partially covered by ice. In spite of their
relevance for primary-productivity estimates arrigo2012massive ;
assmy2017leads , under-ice blooms are still not well understood.
Notwithstanding its idealized character, the model considered here can then be
informative about the basic mechanisms controlling the interplay between
transport and biological growth processes in polar, light-limited
environments, in the prespective of properly representing fluid motions in
more complex numerical models.

Aiming to explore the simplest possible settings, we considered obstacles that
fully absorb light and entirely cover either the downwelling region (LBD case)
or the upwelling one (LBU case). We found that the main role of advective
transport is to hinder phytoplankton growth. In both the LBD and LBU cases,
the per-capita growth rate decreases with increasing flow intensity, and for
sufficiently large advecting velocity, it can be become negative, meaning that
the population goes extinct. While this behavior is qualitatively the same as
the one observed in the absence of horizontal light-intensity modulation, now
the quantitative effect depends on the position of the obstacles with respect
to the regions of ascending and descending fluid motion. Phytoplankton
survival is more difficult when the transmission of light is blocked over the
upwelling flow region (LBU).

The flow field alone cannot be responsible for the asymmetric response of the
system in the LBD and LBU cases, as it can be understood by examining the
advection-reaction dynamics in a Lagrangian reference frame. We then explored
the role of the sinking speed (due to the mass-density difference between
phytoplankton organisms and water). Our analysis shows that, while typically
much smaller than the flow intensity, the sinking speed is at the origin of
the differences observed in the per-capita growth rates. We then provided a
qualitative explanation of this phenomenon, which essentially depends on how
the sinking speed combines with the vertical advecting velocity. In
downwelling regions they sum, while in upwelling regions the sinking speed is
subtracted from the flow. As a consequence, the residence time in the well-lit
upper layer is increased in the upwelling region, relevant for the LBD case,
and decreased in the downwelling region, relevant for the LBU case.

Finally, we have shown that the combination of velocities in the vertical
discussed above is delicate, and can be strongly affected by lateral
diffusion. Indeed, in the present system, vertical and horizontal motions
cannot be easily decoupled. Diffusion spreads the population on the
horizontal, bringing it from regions that are favorable for reproduction to
deadly ones. This effect is more marked in the LBD case, where the population
is more concentrated at the surface (in the upwelling region, where sinking
partially counteracts the flow). In the LBU case, instead, the dynamics are
less sensitive to the intensity of diffusion, because in the (downwelling)
illuminated region the population is already transported (by both the flow and
sinking) to deeper waters where reproduction is also strongly reduced.

To conclude, this study provided evidence that the onset of convection may be
harmful to under-ice blooms, confirming previous findings lowry2018under . The
most notable result, however, concerns the role of the phytoplankton sinking
speed. While small, the latter has a measurable impact on the growth dynamics
of the population and can even be critical for its survival, in partially ice
covered waters undergoing convective motions. This fact may have ecological
relevance, as different phytoplankton species have different densities and,
hence, different settling velocity. The system dynamics revealed highly
complex, also in such an idealized setting, due to the very subtle interplay
bewteen the different transport mechanisms (fluid advection, sinking,
diffusion) and biological growth, which limits the possibility to obtain an
analytical understanding. An interesting perspective would be to carry out
numerical simulations in a more realistic configuration for the fluid, by
considering the dynamics of turbulent convection in the presence of a melting
interface with the ice obstacle esfahani2018basal .

## Acknowledgements

We are grateful to M. Cencini for interesting discussions about the role of
diffusion.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * (1) K. Mann, J. Lazier, Dynamics of marine ecosystems: biological-physical interactions in the oceans (Wiley, New York, 2005)
  * (2) T. Kiørboe, A mechanistic approach to plankton ecology (Princeton University Press, Princeton, 2018)
  * (3) J.C. Prairie, K.R. Sutherland, K.J. Nickols, A.M. Kaltenberg, Limnol. Oceanogr. 2, 121–145 (2012)
  * (4) J.R. Taylor, R. Ferrari, Limnol. Oceanogr. 56, 2293–2307 (2011)
  * (5) C. Lindemann, A. Visser, P. Mariani, Patrizio, J. R. Soc. Interface 14, 20170453 (2017)
  * (6) V.B. Tergolina, E. Calzavarini, G. Mompean, S. Berti, Phys. Rev. E 104, 065106 (2021)
  * (7) K.E. Lowry, R.S. Pickart, V. Selz, M.M. Mills, A. Pacini, K.M. Lewis, H.L. Joy-Warren, C. Nobre, G.L. van Dijken, P.-L. Grondin, J. Ferland, K.R. Arrigo , J. Geophys. Res. 123, 90–109 (2018)
  * (8) M. Ardyna, C.J. Mundy, N. Mayot, L.C. Matthes, L. Oziel, C. Horvat, E. Leu, P. Assmy, V. Hill, P. Matrai, M. Gale, I.A. Melnikov, K.R. Arrigo, Front. Mar. Sci. 7, 608032 (2020)
  * (9) K.R. Arrigo, D.K. Perovich, R.S. Pickart, Z.W. Brown, G.L. van Dijken, K.E. Lowry, M.M. Mills, M.A. Palmer, W.M. Balch, F. Bahr, N.R. Bates, C. Benitez-Nelson, B. Bowler, E. Brownlee, J.K. Ehn, K.E. Frey, R. Garley, S.R. Laney, L. Lubelczyk, J. Mathis, A. Matsuoka, B.G. Mitchell, G.W.K. Moore, E. Ortega-Retuerta, S. Pal, C.M. Polashenski, R.A. Reynolds, B. Schieber, H.M. Sosik, M. Stephens, J.H. Swift, Science 336, 1408 (2012)
  * (10) K.R. Arrigo, G.L. van Dijken, Progr. Oceanogr. 136, 60–70 (2015)
  * (11) P. Assmy, M. Fernández-Méndez, P. Duarte, A. Meyer, A. Randelhoff, C.J. Mundy, L.M. Olsen, H.M. Kauko, A. Bailey, M. Chierici, L. Cohen, A.P. Doulgeris, J.K. Ehn, A. Fransson, S. Gerland, H. Hop, S.R. Hudson, N. Hughes, P. Itkin, G. Johnsen, J.A. King, B.P. Koch, Z. Koenig, S. Kwasniewski, S.R. Laney, M. Nicolaus, A.K. Pavlov, C.M. Polashenski, C. Provost, A. Rösel, M. Sandbu, G. Spreen, L.H. Smedsrud, A. Sundfjord, T. Taskjelle, A. Tatarek, J. Wiktor, P.M. Wagner, A. Wold, H. Steen, M.A. Granskog, Sci. Rep. 7, 1–9 (2017)
  * (12) E. Boles, C. Provost, V. Garçon, C. Bertosio, M. Athanase, Z. Koenig, N. Sennéchael, J. Geophys. Res. 125, e2019JC015608 (2020)
  * (13) H.U. Sverdrup, HU, J. Cons. Int. Explor. Mer 18, 287–295 (1953)
  * (14) R.W. Lindsay, J. Zhang, J. Clim. 18, 4879–4894 (2005)
  * (15) M.C. Serreze, M.M. Holland, J. Stroeve, Julienne, Science 315, 1533–1536 (2007)
  * (16) R. Kwok, G. Spreen, S. Pang, J. Geophys. Res. 118, 2408–2425 (2013)
  * (17) N. Shigesada, A. Okubo, J. Math. Biol. 12, 311–326 (1981)
  * (18) J. Huisman, F.J. Weissing, Ecology 75, 507–520 (1994)
  * (19) J. Huisman, M. nd Arrayás, U. Ebert, B. Sommeijer, Am. Nat. 159, 245–254 (2002)
  * (20) J. Huisman, B. Sommeijer, J. Sea Res. 48, 83–96 (2002)
  * (21) S.A. Thorpe, The turbulent ocean (Cambridge University Press, New York, 2005)
  * (22) A.D.D. Craik, S. Leibovich, J. Fluid Mech. 73, 401–426 (1976)
  * (23) M.A. Bees, Appl. Sci. Res. 59, 141–158 (1997)
  * (24) S.A. Thorpe, Ann. Rev. Fluid Mech. 36, 55–79 (2004)
  * (25) G.K. Vallis, Atmospheric and oceanic fluid dynamics (Cambridge University Press, New York, 2017)
  * (26) J. Huisman, P. van Oostveen, F.J. Weissing, Limnol. Oceanogr. 44, 1781–1787 (1999)
  * (27) T.H. Solomon, J.P. Gollub, Phys. Rev. A 38, 6280 (1988)
  * (28) G. Lacorata, A. Vulpiani, Phys. Rev. E 95, 043106 (2017)
  * (29) G. La Forgia, D. Cavaliere, S. Espa, F. Falcini, G. Lacorata, Sci. Rep., 12, 1–12 (2022)
  * (30) M. Abel, A. Celani, D. Vergni, A. Vulpiani, Phys. Rev. E 64, 046307 (2001)
  * (31) S. Berti, D. Vergni, F. Visconti, A. Vulpiani, Phys. Rev. E 72, 036302 (2005)
  * (32) M. Nicolaus, C. Katlein, J. Maslanik, S. Hendricks, Geophys. Res. Lett. 39, L24501 (2012)
  * (33) H. Stommel, J. Mar. Res. 8, 24–29 (1949)
  * (34) A.B. Ryabov, B. Blasius, Math. Model. Nat. Phenom. 3, 42–86 (2008)
  * (35) B. Rabbanipour Esfahani, S.C. Hirata, S. Berti, E. Calzavarini, Phys. Rev. Fluids 3, 053501 (2018)