Are dependent variables necessarily functions of one another? The problem
Suppose you have two variables $X_1,X_2$ so that $X_1\not\perp\!\!\!\!\! \perp X_2$.
Do we necessarily have that a functional relationship exists between them?
I am assuming random variables are the usual real valued measurable functions, and if that helps one can assume every variable has a density and the conditional densities exist.
I am interested in potential transforms $f\colon\mathbb{R}^k\to\mathbb{R}$ which act on the realisations of one variable (say $X_2$) to (potentially nondeterministically) give the other (say $X_1$).
I am ruling out the cases of confounding and conditioning on colliders.
Formally, do we have that dependence necessarily implies $\exists f,Z\;$ so that $X_1 = f(X_2,Z)$ and either $Z\perp\!\!\! \perp X_2$ or $Z\perp\!\!\! \perp X_1$ .
$Z$ should be interpreted as a disturbance term responsible for the non-deterministic nature of the mapping.
Obviously, I want the function $f$ to be so that $f(x,\cdot)\neq f(x',\cdot)$ for some $x,x'\in \text{Supp}(X_2)$. Intuitively I'm interested in describing $X_1$ as the result of a non-deterministic mapping of $X_2$.
Further, I assume for now that a variable that has the same marginal distribution as $X_1$ and the same joint distribution with $X_2$ can be identified with $X_1$. The idea is "if it behaves the same, why make a distinction?".
The motivation
I essentially would like some tools to justify working with models of the form $X_1 = f(X_2) + Z$ or in a more general case $X_1 = f(X_2,Z)$.
Some ideas
One idea I can think of is the Renyi independence criterion: if $X_1,X_2$ are dependent then for some continuous functions $f,g$ we have $\text{Cov}_{X_1,X_2}(f(X_1),g(X_2))\neq 0$.
This implies that for a region of non-zero probability, $f(X_1) = \alpha g(X_2)+ \beta$.
However this doesn't help me further than this. Playing with inverse cdfs also doesn't really work as far as I'm concerned.
Some examples
These include


*

*additive noise $Z$ on top of some univariate function as $X_1=f(X_2)+Z$. Then $X_1\mid X_2$ is simply the distribution of $Z$ with an altered location $f(x_2)$

*Multiplicative noise, for example @Dave 's answer: $X_1 = Z\cdot f(X_2)$. If $Z$ takes values $\{-1,1\}$ you do observe two curves, but they are described as a unique functional relationship

 A: This question is interesting because it concerns the general problem of regression in the sense of characterizing (or estimating based on data) the conditional distribution of one variable, $X_1,$ based on values of another variable $X_2.$  Their lack of statistical independence implies there is something to be gained from this.

We may take our cues from the theory of Generalized Linear Models (GLMs), which suppose the conditional distribution of $X_1$ is a member of a finitely parameterized family of distributions and the parameters depend first on a linear function of $X_2$ and, usually, by a nonlinear transformation of that result (aka the "link function").  In the present question no such parametric assumption is made.  We should therefore work with the conditional distribution of $X_1$ directly.
Suppose, then, that $X_2=x_2$ is a given value and that we can make sense of the conditional distribution of $X_1$ given $X_2=x_2.$  Let this distribution function be $F_{x_2}.$  (For the moment, let's drop the $x_2$ subscript for brevity.)  By definition, for any number $x_1,$
$$\Pr(X_1 \le x_1\mid X_2=x_2) = F(x_1).$$
The axioms of probability imply $F$ has a right inverse $F^{-1}$ characterized by
$$F\left(F^{-1}(p)\right) = p$$
for all $0\lt p \lt 1.$
Let $Z$ be a random variable with a uniform distribution on the interval $(0,1).$  This means that for all $0\le p\le 1,$ $\Pr(Z \le p) = p$.  Note that for any number $x,$
$$\Pr\left(F^{-1}(Z) \le x\right) = \Pr\left(F(F^{-1}(Z)) \le F(x)\right) = \Pr\left(Z \le F(x)\right) = F(x).$$
This shows that the distribution function of the random variable $F^{-1}(Z)$ is $F.$
Suppose it is possible to find such a uniform $Z$ that is independent of $X_2.$ (One can always do this by enlarging the original probability space if necessary.)  In this case the foregoing result holds for every possible value $x_2.$  Consequently, if we define the function $f:\mathbb{R}^2\to \mathbb{R}$ by

$$f(x_2, z) = F_{x_2}^{-1}(z),$$

then $(f(X_2,Z), X_2)$ is a bivariate random variable for which

*

*Its second component, $X_2,$ is identical to the second component of the original variable $(X_1,X_2);$ and


*The conditional distribution of its first component equals the conditional distribution of $X_1$ given $X_2.$
Consequently,

The random variables $(X_1,X_2)$ and $(f(X_2,Z), X_2)$ are equal in distribution.

A fortiori, $f(X_2,Z)$ and $X_1$ are equal in distribution.

Examples
As the first example, consider the intriguing illustration presented in the reply by Dave at https://stats.stackexchange.com/a/443150/919/.  It displays a $42$-element dataset consisting of a sequence of ordered pairs
$$(x_1,x_2) = (i/10, \pm i/10)$$
as $i$ ranges from $-20$ to $20$ inclusive.  (The value $(0,0)$ is included twice, even though that is not apparent from the scatterplot.)  We may conceive of this in terms of its empirical distribution, which assigns the probability $1/42$ to each element of the dataset, thereby defining a bivariate random variable $(X_1,X_2).$
In this case we will need to define $F_{x_2}$ for each value $x_2$ that appears as a second coordinate in the dataset (it doesn't matter how we define $F_{x_2}$ elsewhere, because those values have no probability). With this in mind, compute
$$F_{i}(x) = \left\{\matrix{0 & x \lt -|i| \\ 1/2 & -|i| \le x \lt |i| \\ 1 & \text{otherwise}}\right.$$
for any number $i.$  One right inverse of $F_{i}$ maps the interval $[0,1/2]$ to  $-|i|$ and the interval $(1/2,1]$ to $|i|.$  You can see this implemented in the code below.
The claim is that

the distribution of $(X_1, X_2)$ is the same as that of $(F_{X_2}^{-1}(Z), X_2).$

Rather than do a direct calculation, let's have the computer simulate the latter using R, to be compared to the scatterplot of $(X_1, X_2)$ in Dave's reply:
set.seed(17)
n <- 1e5                                                  # Simulation size
f <- function(i, p) ifelse (p <= 1/2, -1, 1) * abs(i)     # Inverse conditional distribution
X.2 <- sample((-20):20 / 10, n, replace=TRUE)             # Sample of X[2]
Z <- runif(n)                                             # Independent sample of Z
X.1 <- f(X.2, Z)                                          # Construction of X[1]

This scatterplot confirms it works correctly:

Another example suggested in comments is a uniform distribution on a circle:

Generating these data required only two changes to the code: $(1+X_2)/2$ has a Beta$(1/2,1/2)$ distribution and $F^{-1}$ now reflects the equation of the circle:
X.2 <- rbeta(n, 1/2, 1/2) * 2 - 1
F.inv <- function(p, i) ifelse (p <= 1/2, -1, 1) * sqrt(1 - i^2)

Finally, the more difficult situations occur when the scatterplot cannot be conceived of as exhibiting either of the $X_i$ as a function of the other, as in this uniform distribution on a square:

Here, the marginal distribution of $X_2$ is a mixture of a uniform distribution on the interval $[-1,1]$ and equal point masses at $\pm 1,$ which can be simulated via
X.2 <- pmax(-1, pmin(1, runif(n, -2, 2)))

The function $F^{-1}_{i}$ depends on whether its argument is in the interval $(-1,1)$ or equal to $\pm1;$ it can be implemented as
F.inv <- function(p, i) ifelse (abs(i)==1, 2*p - 1, ifelse (p <= 1/2, -1, 1))

A: I think this is a good place to use a picture. 
x0 <- x1 <- seq(-2,2,0.1)
y0 <- x0
y1 <- -x1
x <- c(x0,x1)
y <- c(y0,y1)
plot(x,y)


I think the lack of independence is evident. However, the relationship cannot be functional, as a large value of $X1$, say $X1=2$, results in either $X2=2$ or $X2=-2$, and there's the usual business about a function having a unique output for each input. 
(I considered that we could call it functional if we see the output as the vector $[-2,2]$ or the set $\{-2,2\}$, but I don't think those work. The relationship is to select one of those values, not a vector of possible values or the entire set of possible values.) 
A: In common factor models, observed variables dependent on the same common factor covary and yet are conditionally independent given the common factor upon which they are jointly dependent. In the literature on graphical causal models, this structure--two variables jointly dependent on a third variable but with no relation between themselves--identifies the third variable as a "confounder," because it produces covariance between variables without a causal relation (see Felix Elwert's very nice introduction).
Alternatively, the two variables may be causes of a third variable.  This marks the third variable as a "collider." Here, conditioning on the third variable (which is a mistake) creates a dependency between the first two, even though there is no functional relation between them.
