Recursive for loop

What I think you are trying to achieve can be done with recursion:

def recursive(n, p=2, prev=[]):
    productList = []
    for k in range(10**(n-1),10**n):
        if p == 1:
            productList.append([k, k])
        else:
            for product in recursive(n, p=p-1):
                productList.append([k] + product[:-1] + [k*product[-1]])
    return productList

do you need all products to be unique, as this will give (2*4=8) also as (4*2=8)?

The batteries for this in SymPy are already included:

from sympy.utilities.iterables import subsets
from sympy.core.mul import prod
numbers = [1, 2, 3]  # or whatever factors you want
p = 2  # how many factors you want
for n in subsets(numbers, p, repetition=True):
    print(n + (prod(n),))

This outputs

(1, 1, 1)
(1, 2, 2)
(1, 3, 3)
(2, 2, 4)
(2, 3, 6)
(3, 3, 9)

This is different from cartes (or product) in that you won't get permutations of the factors, e.g. the latter would also give (2, 1, 2), a permutation of (1, 2, 2). You can read more about subsets using help(subsets) in SymPy and you can read the source if you want to.

The itertools.permutations method is perfect for what you are trying to achieve.

First generate a list of numbers/digits that can be used.

digits = list(range(10 ** (n - 1), 10 ** n))

Then multiply that list by p in order to get p copies of each number.

digits *= p

Finally, use itertools to generate a list of all the permutations (or combinations) possible with a length of p.

for values in itertools.permutations(digits, p):
    output.append([values, np.prod(values)])

The finished function is shown below:

import itertools
import numpy as np

def product(n, p=2):
    output = []

    # Generate a list containing every possible digit.
    digits = list(range(10 ** (n - 1), 10 ** n))

    # Multiply the list by <p> so that we have <p> copies of each digit.
    digits *= p

    # Iterate over each possible permutation of the digits.
    for values in itertools.permutations(digits, p):
        output.append([values, np.prod(values)])

    return output

Alternatively, using numpy and sympy you can make the function much faster (albeit at the cost of having the output be an array rather than list).
Here is the code for a numpy + sympy implementation that is faster (citation needed).

import numpy as np
from sympy.utilities.iterables import multiset_permutations as perm

def product2(n, p=2):
    digits = np.arange(10 ** (n - 1), 10 ** n)
    digits = np.repeat(digits, p)

    permutations = list(perm(digits, size=p))
    output = np.zeros((len(permutations), n + 1))

    for i in range(len(permutations)):
        for j in range(n):
            output[i][j] = permutations[i][j]
        output[i][-1] = np.prod(permutations[i])

    return output