This is an answer to the part of your question concerning the case of dimension $1$. I'll omit any details in higher dimension; even if you're only vaguely familiar with elliptic curves you will see the bigger picture.

1. A Pell conic is an affine curve of the form $C_N: Q_0(X,Y) = 1$, where $Q_0$ is the principal binary form with discriminant $N$ (if $N = 4m$ with $m \equiv 3 \bmod 4$, then $Q(X,Y) = X^2 - mY^2$). For each prime $p$ not dividing $N$, the curve has a smooth reduction modulo $p$, and the number of points is $p - a_p$ for $a_p = (N/p)$ (the Kronecker symbol).

2. We call $C_N$ modular if there exists a modulus $m$ such that $a_p$ only depends on the residue class of $p$ modulo $m$. It can be shown that $m = N$ always works, and that $(N/p) = \chi(p)$ for a primitive Dirichlet character $\chi$ with conductor $N$.

3. We attach a "Pell form" $$ f_N(q) = \sum_{n=1}^\infty a_n q^n $$ to $C_N$. If the Pell conic is modular, then $f_N$ is a rational function of $q$ (the numerators are basically what are known as Fekete polynomials, and a partial fraction decomposition naturally will lead you to Gauss sums) and can be extended to a meromorphic function on the whole complex plane with poles at the primitive $N$-th roots of unity and satisfying a functional equation $$f_N\Big( \frac1q \Big) = - \chi(-1) f_N(q). $$ Observe that if $\chi(-1) = 1$, then $f_N(1) = 0$. This is closely related to the global solvability of the Pell equation in integers. If $\chi(-1) = -1$, on the other hand, then $f_N(1) = 2h/w$, where $h$ is the class number and $w$ the number of roots of unity in the complex quadratic number field with discriminant $-N$.

4. For a prime $p \nmid N$, define the Hecke operator $T_p$ by $$ f|_{T_p}(z) = \frac1p \sum_{a=0}^{p-1} f\Big(\frac{z+a}p \Big). $$ Then the Pell forms $f_N$ are simultaneous eigenforms of the $T_p$ with eigenvalues $a_p$.

There are other eigenfunctions of $T_p$: the cotangent function, the Bernoulli polynomials, Hurwitz's zeta function etc. This is closely related to what Kubert and Lang have called distributions (cf. etc. Washington's cyclotomic fields).

5. The analog of modular curves are the cyclotomic units. You can parametrize Pell conics analytically by trigonometric functions and "arithmetically" by taking norms from cyclotomic units.

I am currently trying to write this (and quite a bit more) up properly. When I'm done, I'll put a link here. For the time being, Darmon has a paper or two on this topic, and Ash's and Gross's "Fearless symmetry" also goes in this direction. Googling "Pell conics" will lead you to a couple of older "preprints".

This is an answer to the part of your question concerning the case of dimension $1$. I'll omit any details in higher dimension; even if you're only vaguely familiar with elliptic curves you will see the bigger picture.

1. A Pell conic is an affine curve of the form $C_N: Q_0(X,Y) = 1$, where $Q_0$ is the principal binary form with discriminant $N$ (if $N = 4m$ with $m \equiv 3 \bmod 4$, then $Q(X,Y) = X^2 - mY^2$). For each prime $p$ not dividing $N$, the curve has a smooth reduction modulo $p$, and the number of points is $p - a_p$ for $a_p = (N/p)$ (the Kronecker symbol).

2. We call $C_N$ modular if there exists a modulus $m$ such that $a_p$ only depends on the residue class of $p$ modulo $m$. It can be shown that $m = N$ always works, and that $(N/p) = \chi(p)$ for a primitive Dirichlet character $\chi$ with conductor $N$.

3. We attach a "Pell form" $$ f_N(q) = \sum_{n=1}^\infty a_n q^n $$ to $C_N$. If the Pell conic is modular, then $f_N$ is a rational function of $q$ (the numerators are basically what are known as Fekete polynomials, and a partial fraction decomposition naturally will lead you to Gauss sums) and can be extended to a meromorphic function on the whole complex plane with poles at the primitive $N$-th roots of unity and satisfying a functional equation $$f_N\Big( \frac1q \Big) = - \chi(-1) f_N(q). $$ Observe that if $\chi(-1) = 1$, then $f_N(1) = 0$. This is closely related to the global solvability of the Pell equation in integers. If $\chi(-1) = -1$, on the other hand, then $f_N(1) = 2h/w$, where $h$ is the class number and $w$ the number of roots of unity in the complex quadratic number field with discriminant $-N$.

4. For a prime $p \nmid N$, define the Hecke operator $T_p$ by $$ f|_{T_p}(z) = \frac1p \sum_{a=0}^{p-1} f\Big(\frac{z+a}p \Big). $$ Then the Pell forms $f_N$ are simultaneous eigenforms of the $T_p$ with eigenvalues $a_p$.

There are other eigenfunctions of $T_p$: the cotangent function, the Bernoulli polynomials, Hurwitz's zeta function etc. This is closely related to what Kubert and Lang have called distributions (cf. etc. Washington's cyclotomic fields).

5. The modular analog of Pell conics are the cyclotomic units. You can parametrize Pell conics analytically by trigonometric functions and "arithmetically" by taking norms from cyclotomic units.

Edit. Here's the link to the preprint.